sc_neurocore_engine/neurons/biophysical/
mihalas_niebur.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 — Mihalas-Niebur Neuron Model
8
9//! Mihalas-Niebur generalized integrate-and-fire dynamics.
10
11/// Mihalas-Niebur 2009 — generalised IF capturing 20 spike patterns.
12#[derive(Clone, Debug)]
13pub struct MihalasNieburNeuron {
14    pub v: f64,
15    pub theta: f64,
16    pub i1: f64,
17    pub i2: f64,
18    pub v_rest: f64,
19    pub v_reset: f64,
20    pub theta_reset: f64,
21    pub theta_inf: f64,
22    pub tau_v: f64,
23    pub tau_theta: f64,
24    pub tau_1: f64,
25    pub tau_2: f64,
26    pub a: f64,
27    pub b: f64,
28    pub r1: f64,
29    pub r2: f64,
30    pub dt: f64,
31}
32
33impl MihalasNieburNeuron {
34    pub fn new() -> Self {
35        Self {
36            v: 0.0,
37            theta: 1.0,
38            i1: 0.0,
39            i2: 0.0,
40            v_rest: 0.0,
41            v_reset: 0.0,
42            theta_reset: 1.0,
43            theta_inf: 1.0,
44            tau_v: 10.0,
45            tau_theta: 100.0,
46            tau_1: 10.0,
47            tau_2: 200.0,
48            a: 0.0,
49            b: 0.0,
50            r1: 0.0,
51            r2: 0.0,
52            dt: 1.0,
53        }
54    }
55
56    fn finite_values(values: &[f64]) -> bool {
57        values.iter().all(|value| value.is_finite())
58    }
59
60    fn valid_runtime(&self) -> bool {
61        Self::finite_values(&[
62            self.v,
63            self.theta,
64            self.i1,
65            self.i2,
66            self.v_rest,
67            self.v_reset,
68            self.theta_reset,
69            self.theta_inf,
70            self.tau_v,
71            self.tau_theta,
72            self.tau_1,
73            self.tau_2,
74            self.a,
75            self.b,
76            self.r1,
77            self.r2,
78            self.dt,
79        ]) && self.tau_v > 0.0
80            && self.tau_theta > 0.0
81            && self.tau_1 > 0.0
82            && self.tau_2 > 0.0
83            && self.dt > 0.0
84    }
85
86    fn derivatives(&self, v: f64, theta: f64, i1: f64, i2: f64, current: f64) -> [f64; 4] {
87        [
88            (-(v - self.v_rest) + i1 + i2 + current) / self.tau_v,
89            (self.theta_inf - theta + self.a * (v - self.v_rest)) / self.tau_theta,
90            -i1 / self.tau_1,
91            -i2 / self.tau_2,
92        ]
93    }
94
95    fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
96        [
97            state[0] + scale * slope[0],
98            state[1] + scale * slope[1],
99            state[2] + scale * slope[2],
100            state[3] + scale * slope[3],
101        ]
102    }
103
104    fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
105        let state = [self.v, self.theta, self.i1, self.i2];
106        let half_dt = 0.5 * self.dt;
107        let k1 = self.derivatives(state[0], state[1], state[2], state[3], current);
108        let s2 = Self::add_scaled(state, k1, half_dt);
109        let k2 = self.derivatives(s2[0], s2[1], s2[2], s2[3], current);
110        let s3 = Self::add_scaled(state, k2, half_dt);
111        let k3 = self.derivatives(s3[0], s3[1], s3[2], s3[3], current);
112        let s4 = Self::add_scaled(state, k3, self.dt);
113        let k4 = self.derivatives(s4[0], s4[1], s4[2], s4[3], current);
114        let candidate = [
115            state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
116            state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
117            state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
118            state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
119        ];
120        if Self::finite_values(&candidate) {
121            Some(candidate)
122        } else {
123            None
124        }
125    }
126
127    pub fn step(&mut self, current: f64) -> i32 {
128        if !current.is_finite() || !self.valid_runtime() {
129            return 0;
130        }
131        let Some(candidate) = self.rk4_candidate(current) else {
132            return 0;
133        };
134        self.v = candidate[0];
135        self.theta = candidate[1];
136        self.i1 = candidate[2];
137        self.i2 = candidate[3];
138        if self.v >= self.theta {
139            self.v = self.v_reset + self.b * (self.v - self.v_rest);
140            self.theta = self.theta_reset.max(self.theta);
141            self.i1 += self.r1;
142            self.i2 += self.r2;
143            1
144        } else {
145            0
146        }
147    }
148
149    /// Run `n_steps` of the candidate-first RK4 recurrence under a constant
150    /// `current`, recording the membrane voltage after every step.
151    ///
152    /// Reuses [`step`] verbatim so the compiled inner loop is bit-identical to
153    /// the per-step path; returns the voltage trace and the total spike count.
154    pub fn simulate(&mut self, n_steps: usize, current: f64) -> (Vec<f64>, i64) {
155        let mut trace = Vec::with_capacity(n_steps);
156        let mut spikes: i64 = 0;
157        for _ in 0..n_steps {
158            spikes += i64::from(self.step(current));
159            trace.push(self.v);
160        }
161        (trace, spikes)
162    }
163
164    pub fn reset(&mut self) {
165        self.v = self.v_rest;
166        self.theta = self.theta_reset;
167        self.i1 = 0.0;
168        self.i2 = 0.0;
169    }
170}
171impl Default for MihalasNieburNeuron {
172    fn default() -> Self {
173        Self::new()
174    }
175}
176
177#[cfg(test)]
178mod tests {
179    use super::*;
180
181    #[test]
182    fn default_matches_constructor_state() {
183        let default = MihalasNieburNeuron::default();
184        let constructed = MihalasNieburNeuron::new();
185        assert_eq!(default.v, constructed.v);
186    }
187
188    #[test]
189    fn simulate_matches_repeated_steps() {
190        let mut simulated = MihalasNieburNeuron::new();
191        let mut repeated = simulated.clone();
192        let (trace, spikes) = simulated.simulate(8, 5.0);
193        let mut expected_spikes = 0_i64;
194        let expected: Vec<f64> = (0..8)
195            .map(|_| {
196                expected_spikes += i64::from(repeated.step(5.0));
197                repeated.v
198            })
199            .collect();
200        assert_eq!(trace, expected);
201        assert_eq!(spikes, expected_spikes);
202    }
203
204    #[test]
205    fn nonfinite_rk4_candidate_preserves_state() {
206        let mut n = MihalasNieburNeuron::new();
207        n.dt = f64::MAX;
208        let before = (n.v, n.theta, n.i1, n.i2);
209        assert_eq!(n.step(1.0), 0);
210        assert_eq!((n.v, n.theta, n.i1, n.i2), before);
211    }
212
213    #[test]
214    fn mn_fires() {
215        let mut n = MihalasNieburNeuron::new();
216        let t: i32 = (0..100).map(|_| n.step(5.0)).sum();
217        assert!(t > 0);
218    }
219
220    // -- MihalasNiebur --
221    #[test]
222    fn mn_silent_without_input() {
223        let mut n = MihalasNieburNeuron::new();
224        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
225        assert_eq!(t, 0);
226    }
227    #[test]
228    fn mn_reset_clears_state() {
229        let mut n = MihalasNieburNeuron::new();
230        for _ in 0..100 {
231            n.step(5.0);
232        }
233        n.reset();
234        assert!((n.v - n.v_rest).abs() < 1e-10);
235        assert!((n.theta - n.theta_reset).abs() < 1e-10);
236    }
237    #[test]
238    fn mn_rk4_reference_point() {
239        let mut n = MihalasNieburNeuron::new();
240        assert_eq!(n.step(0.5), 0);
241        assert!((n.v - 0.04758125).abs() < 1e-12);
242        assert!((n.theta - 1.0).abs() < 1e-15);
243        assert!((n.i1 - 0.0).abs() < 1e-15);
244        assert!((n.i2 - 0.0).abs() < 1e-15);
245    }
246    #[test]
247    fn mn_spike_reset_uses_b() {
248        let mut n = MihalasNieburNeuron::new();
249        n.v = 0.99;
250        n.b = 0.5;
251        n.r1 = 1.25;
252        n.r2 = -0.25;
253        assert_eq!(n.step(2.0), 1);
254        assert!((n.v - 0.5430570625).abs() < 1e-12);
255        assert!((n.i1 - 1.25).abs() < 1e-15);
256        assert!((n.i2 - (-0.25)).abs() < 1e-15);
257    }
258    #[test]
259    fn mn_invalid_input_preserves_state() {
260        let mut n = MihalasNieburNeuron::new();
261        n.v = 0.2;
262        n.i1 = 0.3;
263        let before = (n.v, n.theta, n.i1, n.i2);
264        assert_eq!(n.step(f64::NAN), 0);
265        assert_eq!((n.v, n.theta, n.i1, n.i2), before);
266    }
267    #[test]
268    fn mn_extreme_bounded() {
269        let mut n = MihalasNieburNeuron::new();
270        for _ in 0..200 {
271            n.step(1e4);
272        }
273        assert!(n.v.is_finite());
274    }
275    #[test]
276    fn mn_adaptive_threshold() {
277        let mut n = MihalasNieburNeuron::new();
278        n.a = 0.1;
279        for _ in 0..100 {
280            n.step(5.0);
281        }
282        // Threshold should have adapted
283        assert!(n.theta.is_finite());
284    }
285    #[test]
286    fn mn_negative_no_crash() {
287        let mut n = MihalasNieburNeuron::new();
288        for _ in 0..200 {
289            n.step(-10.0);
290        }
291        assert!(n.v.is_finite());
292    }
293    #[test]
294    fn mn_nan_no_panic() {
295        let mut n = MihalasNieburNeuron::new();
296        n.step(f64::NAN);
297    }
298}
sc_neurocore_engine/neurons/biophysical/mihalas_niebur.rs

sc_neurocore_engine/neurons/biophysical/
mihalas_niebur.rs
