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

sc_neurocore_engine/neurons/biophysical/
glif.rs
