sc_neurocore_engine/neurons/biophysical/
gif_population.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 — GIF Population Neuron Model
8
9//! Seeded escape-rate generalized integrate-and-fire population dynamics.
10
11use rand::{RngExt, SeedableRng};
12use rand_xoshiro::Xoshiro256PlusPlus;
13
14/// GIF population — escape-rate generalized IF. Mensi et al. 2012.
15#[derive(Clone, Debug)]
16pub struct GIFPopulationNeuron {
17    pub v: f64,
18    pub theta: f64,
19    pub eta: f64,
20    pub tau_m: f64,
21    pub tau_eta: f64,
22    pub delta_v: f64,
23    pub lambda_0: f64,
24    pub eta_increment: f64,
25    pub v_rest: f64,
26    pub v_reset: f64,
27    pub dt: f64,
28    pub seed: u64,
29    rng: Xoshiro256PlusPlus,
30}
31
32impl GIFPopulationNeuron {
33    pub fn new(seed: u64) -> Self {
34        Self {
35            v: -65.0,
36            theta: -50.0,
37            eta: 0.0,
38            tau_m: 20.0,
39            tau_eta: 100.0,
40            delta_v: 2.0,
41            lambda_0: 0.001,
42            eta_increment: 5.0,
43            v_rest: -65.0,
44            v_reset: -65.0,
45            dt: 0.5,
46            seed,
47            rng: Xoshiro256PlusPlus::seed_from_u64(seed),
48        }
49    }
50
51    fn finite_values(values: &[f64]) -> bool {
52        values.iter().all(|value| value.is_finite())
53    }
54
55    fn valid_runtime(&self) -> bool {
56        Self::finite_values(&[
57            self.v,
58            self.theta,
59            self.eta,
60            self.tau_m,
61            self.tau_eta,
62            self.delta_v,
63            self.lambda_0,
64            self.eta_increment,
65            self.v_rest,
66            self.v_reset,
67            self.dt,
68        ]) && self.tau_m > 0.0
69            && self.tau_eta > 0.0
70            && self.delta_v > 0.0
71            && self.lambda_0 >= 0.0
72            && self.dt > 0.0
73    }
74
75    fn advance_subthreshold(&self, current: f64) -> Option<(f64, f64)> {
76        let eta_decay = (-self.dt / self.tau_eta).exp();
77        let membrane_decay = (-self.dt / self.tau_m).exp();
78        let x0 = self.v - self.v_rest - current;
79        let eta_new = self.eta * eta_decay;
80        let x_new = if (self.tau_m - self.tau_eta).abs() <= 1e-12 {
81            membrane_decay * (x0 - self.eta * self.dt / self.tau_m)
82        } else {
83            let coupling = self.tau_eta / (self.tau_eta - self.tau_m);
84            x0 * membrane_decay - self.eta * coupling * (eta_decay - membrane_decay)
85        };
86        let v_new = self.v_rest + current + x_new;
87        if Self::finite_values(&[v_new, eta_new]) {
88            Some((v_new, eta_new))
89        } else {
90            None
91        }
92    }
93
94    fn spike_probability(&self, voltage: f64) -> f64 {
95        if self.lambda_0 == 0.0 {
96            return 0.0;
97        }
98        let exponent = ((voltage - self.theta) / self.delta_v).clamp(-745.0, 20.0);
99        let hazard = self.lambda_0 * exponent.exp();
100        (1.0 - (-hazard * self.dt).exp()).clamp(0.0, 1.0)
101    }
102
103    pub fn step(&mut self, current: f64) -> i32 {
104        if !current.is_finite() || !self.valid_runtime() {
105            return 0;
106        }
107        let Some((v_candidate, eta_candidate)) = self.advance_subthreshold(current) else {
108            return 0;
109        };
110        self.v = v_candidate;
111        self.eta = eta_candidate;
112        if self.rng.random::<f64>() < self.spike_probability(self.v) {
113            self.v = self.v_reset;
114            self.eta += self.eta_increment;
115            1
116        } else {
117            0
118        }
119    }
120
121    pub fn reset(&mut self) {
122        self.v = self.v_rest;
123        self.eta = 0.0;
124        self.rng = Xoshiro256PlusPlus::seed_from_u64(self.seed);
125    }
126}
127
128#[cfg(test)]
129mod tests {
130    use super::*;
131
132    #[test]
133    fn equal_time_constants_use_limit_solution() {
134        let mut n = GIFPopulationNeuron::new(42);
135        n.tau_eta = n.tau_m;
136        assert_eq!(n.step(0.0), 0);
137        assert!(n.v.is_finite());
138        assert!(n.eta.is_finite());
139    }
140
141    #[test]
142    fn zero_escape_rate_remains_subthreshold() {
143        let mut n = GIFPopulationNeuron::new(42);
144        n.lambda_0 = 0.0;
145        assert_eq!(n.step(30.0), 0);
146        assert_ne!(n.v, n.v_reset);
147    }
148
149    #[test]
150    fn nonfinite_subthreshold_candidate_preserves_state() {
151        let mut n = GIFPopulationNeuron::new(42);
152        n.v = f64::MAX;
153        n.v_rest = -f64::MAX;
154        let before = (n.v, n.eta);
155        assert_eq!(n.step(f64::MAX), 0);
156        assert_eq!((n.v, n.eta), before);
157    }
158
159    #[test]
160    fn gif_pop_fires() {
161        let mut n = GIFPopulationNeuron::new(42);
162        let t: i32 = (0..1000).map(|_| n.step(30.0)).sum();
163        assert!(t > 0);
164    }
165
166    // -- GIFPopulation --
167    #[test]
168    fn gif_pop_exact_subthreshold_reference_point() {
169        let mut n = GIFPopulationNeuron::new(7);
170        n.v = -68.0;
171        n.eta = 0.4;
172        assert_eq!(n.step(4.0), 0);
173        assert!((n.v - (-67.8370206677805)).abs() < 1e-12);
174        assert!((n.eta - 0.398004991677073).abs() < 1e-15);
175    }
176    #[test]
177    fn gif_pop_forced_spike_adds_decayed_adaptation() {
178        let mut n = GIFPopulationNeuron::new(42);
179        n.v = -51.0;
180        n.eta = 0.3;
181        n.theta = -90.0;
182        n.lambda_0 = 1.0e9;
183        assert_eq!(n.step(0.0), 1);
184        assert!((n.v - n.v_reset).abs() < 1e-12);
185        assert!((n.eta - 5.298503743757805).abs() < 1e-15);
186    }
187    #[test]
188    fn gif_pop_invalid_input_preserves_state() {
189        let mut n = GIFPopulationNeuron::new(42);
190        n.v = -62.0;
191        n.eta = 0.75;
192        let before = (n.v, n.eta);
193        assert_eq!(n.step(f64::NAN), 0);
194        n.tau_m = 0.0;
195        assert_eq!(n.step(1.0), 0);
196        assert_eq!((n.v, n.eta), before);
197    }
198    #[test]
199    fn gif_pop_seeded_reset_replays() {
200        let mut n = GIFPopulationNeuron::new(123);
201        n.theta = -90.0;
202        n.lambda_0 = 1.0e9;
203        let first: Vec<i32> = (0..3).map(|_| n.step(0.0)).collect();
204        n.reset();
205        let replay: Vec<i32> = (0..3).map(|_| n.step(0.0)).collect();
206        assert_eq!(first, replay);
207        assert!(n.eta > n.eta_increment);
208    }
209    #[test]
210    fn gif_pop_negative_drive_remains_finite() {
211        let mut n = GIFPopulationNeuron::new(42);
212        for _ in 0..200 {
213            n.step(-30.0);
214        }
215        assert!(n.v.is_finite());
216        assert!(n.eta.is_finite());
217    }
218}
sc_neurocore_engine/neurons/biophysical/gif_population.rs

sc_neurocore_engine/neurons/biophysical/
gif_population.rs
