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 — source-faithful Mihalas-Niebur generalised IF model
8
9//! Mihalaş-Niebur equations 2.1–2.2 with capacitance-normalised currents.
10
11/// Four-state generalised integrate-and-fire neuron from Mihalaş and Niebur (2009).
12///
13/// Rates are per millisecond, voltages are volts, and currents are volts per
14/// millisecond after division by capacitance. The flow uses fixed-grid RK4 and
15/// sampled threshold detection; the published differential equations and event
16/// reset are otherwise unchanged.
17#[derive(Clone, Debug)]
18pub struct MihalasNieburNeuron {
19    pub v: f64,
20    pub theta: f64,
21    pub i1: f64,
22    pub i2: f64,
23    pub v_rest: f64,
24    pub v_reset: f64,
25    pub theta_reset: f64,
26    pub theta_inf: f64,
27    pub leak_rate: f64,
28    pub threshold_voltage_coupling: f64,
29    pub threshold_decay_rate: f64,
30    pub current_decay_rate_1: f64,
31    pub current_decay_rate_2: f64,
32    pub current_retention_1: f64,
33    pub current_retention_2: f64,
34    pub current_jump_1: f64,
35    pub current_jump_2: f64,
36    pub dt: f64,
37}
38
39impl MihalasNieburNeuron {
40    /// Construct the paper's common Table 1 profile with Figure 1C coupling.
41    pub fn new() -> Self {
42        Self {
43            v: -0.07,
44            theta: -0.05,
45            i1: 0.0,
46            i2: 0.0,
47            v_rest: -0.07,
48            v_reset: -0.07,
49            theta_reset: -0.06,
50            theta_inf: -0.05,
51            leak_rate: 0.05,
52            threshold_voltage_coupling: 0.005,
53            threshold_decay_rate: 0.01,
54            current_decay_rate_1: 0.2,
55            current_decay_rate_2: 0.02,
56            current_retention_1: 0.0,
57            current_retention_2: 1.0,
58            current_jump_1: 0.0,
59            current_jump_2: 0.0,
60            dt: 0.1,
61        }
62    }
63
64    fn finite(values: &[f64]) -> bool {
65        values.iter().all(|value| value.is_finite())
66    }
67
68    fn valid(&self) -> bool {
69        Self::finite(&[
70            self.v,
71            self.theta,
72            self.i1,
73            self.i2,
74            self.v_rest,
75            self.v_reset,
76            self.theta_reset,
77            self.theta_inf,
78            self.leak_rate,
79            self.threshold_voltage_coupling,
80            self.threshold_decay_rate,
81            self.current_decay_rate_1,
82            self.current_decay_rate_2,
83            self.current_retention_1,
84            self.current_retention_2,
85            self.current_jump_1,
86            self.current_jump_2,
87            self.dt,
88        ]) && self.leak_rate > 0.0
89            && self.threshold_decay_rate > 0.0
90            && self.current_decay_rate_1 > 0.0
91            && self.current_decay_rate_2 > 0.0
92            && self.dt > 0.0
93            && self.theta_reset > self.v_reset
94    }
95
96    fn derivatives(&self, state: [f64; 4], current: f64) -> [f64; 4] {
97        [
98            current + state[2] + state[3] - self.leak_rate * (state[0] - self.v_rest),
99            self.threshold_voltage_coupling * (state[0] - self.v_rest)
100                - self.threshold_decay_rate * (state[1] - self.theta_inf),
101            -self.current_decay_rate_1 * state[2],
102            -self.current_decay_rate_2 * state[3],
103        ]
104    }
105
106    fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
107        [
108            state[0] + scale * slope[0],
109            state[1] + scale * slope[1],
110            state[2] + scale * slope[2],
111            state[3] + scale * slope[3],
112        ]
113    }
114
115    fn candidate(&self, current: f64) -> Option<(Self, i32)> {
116        if !self.valid() || !current.is_finite() {
117            return None;
118        }
119        let state = [self.v, self.theta, self.i1, self.i2];
120        let half_dt = 0.5 * self.dt;
121        let k1 = self.derivatives(state, current);
122        let k2 = self.derivatives(Self::add_scaled(state, k1, half_dt), current);
123        let k3 = self.derivatives(Self::add_scaled(state, k2, half_dt), current);
124        let k4 = self.derivatives(Self::add_scaled(state, k3, self.dt), current);
125        let values = [
126            state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
127            state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
128            state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
129            state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
130        ];
131        if !Self::finite(&values) {
132            return None;
133        }
134        let mut next = Self {
135            v: values[0],
136            theta: values[1],
137            i1: values[2],
138            i2: values[3],
139            ..self.clone()
140        };
141        let event = i32::from(next.v >= next.theta);
142        if event == 1 {
143            next.i1 = self.current_retention_1 * next.i1 + self.current_jump_1;
144            next.i2 = self.current_retention_2 * next.i2 + self.current_jump_2;
145            next.v = self.v_reset;
146            next.theta = self.theta_reset.max(next.theta);
147        }
148        next.valid().then_some((next, event))
149    }
150
151    /// Advance one sampled interval, returning `None` for an invalid candidate.
152    pub fn try_step(&mut self, current: f64) -> Option<i32> {
153        let (candidate, event) = self.candidate(current)?;
154        *self = candidate;
155        Some(event)
156    }
157
158    /// Advance one sampled interval and leave state unchanged on invalid input.
159    pub fn step(&mut self, current: f64) -> i32 {
160        self.try_step(current).unwrap_or(0)
161    }
162
163    /// Simulate a constant-current trajectory atomically.
164    pub fn try_simulate(&mut self, n_steps: usize, current: f64) -> Option<(Vec<f64>, i64)> {
165        let mut candidate = self.clone();
166        let mut trace = Vec::with_capacity(n_steps);
167        let mut events = 0_i64;
168        for _ in 0..n_steps {
169            events += i64::from(candidate.try_step(current)?);
170            trace.push(candidate.v);
171        }
172        *self = candidate;
173        Some((trace, events))
174    }
175
176    /// Restore the paper-profile resting state.
177    pub fn reset(&mut self) {
178        self.v = self.v_rest;
179        self.theta = self.theta_inf;
180        self.i1 = 0.0;
181        self.i2 = 0.0;
182    }
183}
184
185impl Default for MihalasNieburNeuron {
186    fn default() -> Self {
187        Self::new()
188    }
189}
190
191#[cfg(test)]
192mod tests {
193    use super::*;
194
195    #[test]
196    fn paper_profile_adapts_under_constant_drive() {
197        let mut neuron = MihalasNieburNeuron::new();
198        let (_, events) = neuron
199            .try_simulate(2500, 0.002)
200            .expect("valid source profile");
201        assert_eq!(events, 13);
202        assert!(neuron.theta > neuron.theta_inf);
203    }
204
205    #[test]
206    fn event_uses_published_reset_map() {
207        let mut neuron = MihalasNieburNeuron::new();
208        neuron.v = -0.049;
209        neuron.i1 = 0.003;
210        neuron.i2 = -0.001;
211        neuron.current_retention_1 = 0.25;
212        neuron.current_retention_2 = 0.5;
213        neuron.current_jump_1 = 0.004;
214        neuron.current_jump_2 = 0.002;
215        assert_eq!(neuron.step(0.02), 1);
216        assert_eq!(neuron.v, neuron.v_reset);
217        assert!(neuron.theta >= neuron.theta_reset);
218        assert!(neuron.i1 > neuron.current_jump_1);
219        assert!(neuron.i2 > 0.0);
220    }
221
222    #[test]
223    fn invalid_batch_is_failure_atomic() {
224        let mut neuron = MihalasNieburNeuron::new();
225        let before = neuron.clone();
226        assert!(neuron.try_simulate(4, f64::NAN).is_none());
227        assert_eq!(neuron.v, before.v);
228        assert_eq!(neuron.theta, before.theta);
229    }
230}
sc_neurocore_engine/neurons/biophysical/mihalas_niebur.rs

sc_neurocore_engine/neurons/biophysical/
mihalas_niebur.rs
