sc_neurocore_engine/neurons/biophysical/
glif.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 — Teeter 2018 GLIF5 source model
8
9//! Five-state GLIF5 dynamics from Teeter et al. (2018), equations 1–8.
10
11/// Teeter GLIF5 membrane, threshold, after-spike-current, and reset state.
12#[derive(Clone, Debug)]
13pub struct GLIFNeuron {
14    pub v: f64,
15    pub theta_spike: f64,
16    pub i_asc1: f64,
17    pub i_asc2: f64,
18    pub theta_voltage: f64,
19    pub refractory_remaining: f64,
20    pub e_l: f64,
21    pub capacitance: f64,
22    pub resistance: f64,
23    pub theta_inf: f64,
24    pub b_spike: f64,
25    pub b_voltage: f64,
26    pub a_voltage: f64,
27    pub k_asc1: f64,
28    pub k_asc2: f64,
29    pub f_v: f64,
30    pub delta_v: f64,
31    pub delta_theta_spike: f64,
32    pub f_asc1: f64,
33    pub f_asc2: f64,
34    pub delta_i_asc1: f64,
35    pub delta_i_asc2: f64,
36    pub refractory_period: f64,
37    pub dt: f64,
38}
39
40impl GLIFNeuron {
41    /// Construct the source-consistent normalized operating profile.
42    pub fn new() -> Self {
43        Self {
44            v: -70.0,
45            theta_spike: 0.0,
46            i_asc1: 0.0,
47            i_asc2: 0.0,
48            theta_voltage: 0.0,
49            refractory_remaining: 0.0,
50            e_l: -70.0,
51            capacitance: 10.0,
52            resistance: 1.0,
53            theta_inf: -50.0,
54            b_spike: 0.01,
55            b_voltage: 0.01,
56            a_voltage: 0.0001,
57            k_asc1: 0.1,
58            k_asc2: 0.005,
59            f_v: 0.0,
60            delta_v: 0.0,
61            delta_theta_spike: 2.0,
62            f_asc1: 1.0,
63            f_asc2: 1.0,
64            delta_i_asc1: 1.0,
65            delta_i_asc2: 0.5,
66            refractory_period: 2.0,
67            dt: 1.0,
68        }
69    }
70
71    fn finite(values: &[f64]) -> bool {
72        values.iter().all(|value| value.is_finite())
73    }
74
75    fn valid(&self) -> bool {
76        Self::finite(&[
77            self.v,
78            self.theta_spike,
79            self.i_asc1,
80            self.i_asc2,
81            self.theta_voltage,
82            self.refractory_remaining,
83            self.e_l,
84            self.capacitance,
85            self.resistance,
86            self.theta_inf,
87            self.b_spike,
88            self.b_voltage,
89            self.a_voltage,
90            self.k_asc1,
91            self.k_asc2,
92            self.f_v,
93            self.delta_v,
94            self.delta_theta_spike,
95            self.f_asc1,
96            self.f_asc2,
97            self.delta_i_asc1,
98            self.delta_i_asc2,
99            self.refractory_period,
100            self.dt,
101        ]) && self.capacitance > 0.0
102            && self.resistance > 0.0
103            && self.b_spike > 0.0
104            && self.b_voltage > 0.0
105            && self.k_asc1 > 0.0
106            && self.k_asc2 > 0.0
107            && self.dt > 0.0
108            && self.refractory_remaining >= 0.0
109            && self.refractory_period >= 0.0
110    }
111
112    fn decay(rate: f64, dt: f64) -> f64 {
113        (-rate * dt).exp()
114    }
115
116    fn exponential_convolution(decay_rate: f64, forcing_rate: f64, dt: f64) -> f64 {
117        let difference = decay_rate - forcing_rate;
118        let scale = 1.0_f64.max(decay_rate.abs()).max(forcing_rate.abs());
119        if difference.abs() <= 1e-12 * scale {
120            dt * (-decay_rate * dt).exp()
121        } else {
122            ((-forcing_rate * dt).exp() - (-decay_rate * dt).exp()) / difference
123        }
124    }
125
126    fn candidate(&self, current: f64) -> Option<(Self, i32)> {
127        if !self.valid() || !current.is_finite() {
128            return None;
129        }
130        if self.refractory_remaining > 0.0 {
131            let mut next = self.clone();
132            next.refractory_remaining = (self.refractory_remaining - self.dt).max(0.0);
133            return Some((next, 0));
134        }
135
136        let total_current = current + self.i_asc1 + self.i_asc2;
137        let membrane_rate = 1.0 / (self.resistance * self.capacitance);
138        let membrane_decay = Self::decay(membrane_rate, self.dt);
139        let equilibrium_offset = self.resistance * total_current;
140        let voltage_offset = self.v - self.e_l;
141        let next_offset =
142            equilibrium_offset + (voltage_offset - equilibrium_offset) * membrane_decay;
143        let next_v = self.e_l + next_offset;
144        let next_theta_spike = self.theta_spike * Self::decay(self.b_spike, self.dt);
145        let next_i_asc1 = self.i_asc1 * Self::decay(self.k_asc1, self.dt);
146        let next_i_asc2 = self.i_asc2 * Self::decay(self.k_asc2, self.dt);
147        let threshold_forcing = equilibrium_offset * (1.0 - Self::decay(self.b_voltage, self.dt))
148            / self.b_voltage
149            + (voltage_offset - equilibrium_offset)
150                * Self::exponential_convolution(self.b_voltage, membrane_rate, self.dt);
151        let next_theta_voltage = self.theta_voltage * Self::decay(self.b_voltage, self.dt)
152            + self.a_voltage * threshold_forcing;
153        let mut next = Self {
154            v: next_v,
155            theta_spike: next_theta_spike,
156            i_asc1: next_i_asc1,
157            i_asc2: next_i_asc2,
158            theta_voltage: next_theta_voltage,
159            refractory_remaining: 0.0,
160            ..self.clone()
161        };
162        if !next.valid() {
163            return None;
164        }
165        if next.v <= self.theta_inf + next.theta_spike + next.theta_voltage {
166            return Some((next, 0));
167        }
168        next.v = self.e_l + self.f_v * (next.v - self.e_l) - self.delta_v;
169        next.theta_spike += self.delta_theta_spike;
170        next.i_asc1 = self.f_asc1 * next.i_asc1 + self.delta_i_asc1;
171        next.i_asc2 = self.f_asc2 * next.i_asc2 + self.delta_i_asc2;
172        next.refractory_remaining = self.refractory_period;
173        next.valid().then_some((next, 1))
174    }
175
176    /// Checked source update; invalid input leaves state unchanged.
177    pub fn try_step(&mut self, current: f64) -> Option<i32> {
178        let (candidate, event) = self.candidate(current)?;
179        *self = candidate;
180        Some(event)
181    }
182
183    /// Network-runner-compatible update with fail-closed invalid-input behavior.
184    pub fn step(&mut self, current: f64) -> i32 {
185        self.try_step(current).unwrap_or(0)
186    }
187
188    /// Run a failure-atomic constant-current batch.
189    pub fn try_simulate(&mut self, n_steps: usize, current: f64) -> Option<(Vec<f64>, i64)> {
190        let mut candidate = self.clone();
191        let mut trace = Vec::with_capacity(n_steps);
192        let mut events = 0_i64;
193        for _ in 0..n_steps {
194            events += i64::from(candidate.try_step(current)?);
195            trace.push(candidate.v);
196        }
197        *self = candidate;
198        Some((trace, events))
199    }
200
201    /// Restore the normalized source-profile state.
202    pub fn reset(&mut self) {
203        self.v = self.e_l;
204        self.theta_spike = 0.0;
205        self.i_asc1 = 0.0;
206        self.i_asc2 = 0.0;
207        self.theta_voltage = 0.0;
208        self.refractory_remaining = 0.0;
209    }
210}
211
212impl Default for GLIFNeuron {
213    fn default() -> Self {
214        Self::new()
215    }
216}
217
218#[cfg(test)]
219mod tests {
220    use super::*;
221
222    #[test]
223    fn default_profile_has_source_five_state_contract() {
224        let neuron = GLIFNeuron::new();
225        assert!(neuron.valid());
226        assert_eq!(
227            neuron.theta_inf + neuron.theta_spike + neuron.theta_voltage,
228            -50.0
229        );
230    }
231
232    #[test]
233    fn reset_rules_and_refractory_state_are_explicit() {
234        let mut neuron = GLIFNeuron::new();
235        neuron.v = -50.0;
236        assert_eq!(neuron.try_step(40.0), Some(1));
237        assert_eq!(neuron.v, -70.0);
238        assert_eq!(neuron.theta_spike, 2.0);
239        assert_eq!(neuron.i_asc1, 1.0);
240        assert_eq!(neuron.i_asc2, 0.5);
241        assert_eq!(neuron.refractory_remaining, 2.0);
242    }
243
244    #[test]
245    fn invalid_batch_is_failure_atomic() {
246        let mut neuron = GLIFNeuron::new();
247        let before = neuron.clone();
248        assert!(neuron.try_simulate(4, f64::NAN).is_none());
249        assert_eq!(neuron.v, before.v);
250        assert_eq!(neuron.theta_spike, before.theta_spike);
251    }
252
253    #[test]
254    fn source_profile_regimes_are_pinned() {
255        for (current, expected) in [(0.0, 0), (22.0, 22), (30.0, 49), (50.0, 80)] {
256            let mut neuron = GLIFNeuron::new();
257            let (_, events) = neuron.try_simulate(1000, current).expect("valid profile");
258            assert_eq!(events, expected, "current={current}");
259        }
260    }
261}
sc_neurocore_engine/neurons/biophysical/glif.rs

sc_neurocore_engine/neurons/biophysical/
glif.rs
