sc_neurocore_engine/neurons/channels/
nmda.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 — Wang 1999 NMDA-autapse pyramidal neuron
8
9/// Wang (1999) pyramidal LIF neuron with two-stage NMDA-autapse kinetics.
10#[derive(Clone, Debug)]
11pub struct NMDANeuron {
12    pub v: f64,
13    pub x_nmda: f64,
14    pub s_nmda: f64,
15    pub ca: f64,
16    pub refractory_remaining: f64,
17    pub c_m: f64,
18    pub g_l: f64,
19    pub v_l: f64,
20    pub g_nmda: f64,
21    pub e_nmda: f64,
22    pub mg_conc: f64,
23    pub alpha_x: f64,
24    pub tau_x: f64,
25    pub alpha_s: f64,
26    pub tau_s: f64,
27    pub kinetic_scale: f64,
28    pub g_ahp: f64,
29    pub v_k: f64,
30    pub alpha_ca: f64,
31    pub tau_ca: f64,
32    pub dt: f64,
33    pub v_threshold: f64,
34    pub v_reset: f64,
35    pub refractory_period: f64,
36}
37
38impl NMDANeuron {
39    pub fn new() -> Self {
40        Self {
41            v: -70.0,
42            x_nmda: 0.0,
43            s_nmda: 0.0,
44            ca: 0.0,
45            refractory_remaining: 0.0,
46            c_m: 0.5,
47            g_l: 0.025,
48            v_l: -70.0,
49            g_nmda: 0.1,
50            e_nmda: 0.0,
51            mg_conc: 1.0,
52            alpha_x: 1.0,
53            tau_x: 2.0,
54            alpha_s: 1.0,
55            tau_s: 80.0,
56            kinetic_scale: 1.0,
57            g_ahp: 0.0,
58            v_k: -85.0,
59            alpha_ca: 0.2,
60            tau_ca: 80.0,
61            dt: 0.05,
62            v_threshold: -52.0,
63            v_reset: -59.0,
64            refractory_period: 2.0,
65        }
66    }
67
68    fn valid(&self) -> bool {
69        let finite = [
70            self.v,
71            self.x_nmda,
72            self.s_nmda,
73            self.ca,
74            self.refractory_remaining,
75            self.c_m,
76            self.g_l,
77            self.v_l,
78            self.g_nmda,
79            self.e_nmda,
80            self.mg_conc,
81            self.alpha_x,
82            self.tau_x,
83            self.alpha_s,
84            self.tau_s,
85            self.kinetic_scale,
86            self.g_ahp,
87            self.v_k,
88            self.alpha_ca,
89            self.tau_ca,
90            self.dt,
91            self.v_threshold,
92            self.v_reset,
93            self.refractory_period,
94        ]
95        .into_iter()
96        .all(f64::is_finite);
97        finite
98            && (-120.0..=80.0).contains(&self.v)
99            && self.x_nmda >= 0.0
100            && (0.0..=1.0).contains(&self.s_nmda)
101            && self.ca >= 0.0
102            && (0.0..=self.refractory_period).contains(&self.refractory_remaining)
103            && (0.01..=10.0).contains(&self.c_m)
104            && (0.0..=1.0).contains(&self.g_l)
105            && (-100.0..=-40.0).contains(&self.v_l)
106            && (0.0..=2.0).contains(&self.g_nmda)
107            && (-10.0..=10.0).contains(&self.e_nmda)
108            && (0.0..=5.0).contains(&self.mg_conc)
109            && (0.0..=10.0).contains(&self.alpha_x)
110            && (0.01..=100.0).contains(&self.tau_x)
111            && (0.0..=10.0).contains(&self.alpha_s)
112            && (1.0..=1000.0).contains(&self.tau_s)
113            && (0.01..=100.0).contains(&self.kinetic_scale)
114            && (0.0..=10.0).contains(&self.g_ahp)
115            && (-120.0..=-40.0).contains(&self.v_k)
116            && (0.0..=10.0).contains(&self.alpha_ca)
117            && (1.0..=1000.0).contains(&self.tau_ca)
118            && self.dt > 0.0
119            && self.dt <= 0.05
120            && (-80.0..=-30.0).contains(&self.v_threshold)
121            && self.v_reset >= -100.0
122            && self.v_reset < self.v_threshold
123            && (0.0..=20.0).contains(&self.refractory_period)
124    }
125
126    fn derivatives(&self, v: f64, x_nmda: f64, s_nmda: f64, ca: f64, current: f64) -> [f64; 4] {
127        let mg_block = 1.0 / (1.0 + self.mg_conc * (-0.062 * v).exp() / 3.57);
128        let i_l = self.g_l * (v - self.v_l);
129        let i_ahp = self.g_ahp * ca * (v - self.v_k);
130        let i_nmda = self.g_nmda * s_nmda * mg_block * (v - self.e_nmda);
131        [
132            (-i_l - i_ahp - i_nmda + current) / self.c_m,
133            self.kinetic_scale * (-x_nmda / self.tau_x),
134            self.kinetic_scale * (self.alpha_s * x_nmda * (1.0 - s_nmda) - s_nmda / self.tau_s),
135            -ca / self.tau_ca,
136        ]
137    }
138
139    fn rk2_candidate(&self, v: f64, current: f64) -> [f64; 4] {
140        let state = [v, self.x_nmda, self.s_nmda, self.ca];
141        let k1 = self.derivatives(state[0], state[1], state[2], state[3], current);
142        let half_dt = 0.5 * self.dt;
143        let midpoint = [
144            state[0] + half_dt * k1[0],
145            state[1] + half_dt * k1[1],
146            state[2] + half_dt * k1[2],
147            state[3] + half_dt * k1[3],
148        ];
149        let k2 = self.derivatives(midpoint[0], midpoint[1], midpoint[2], midpoint[3], current);
150        [
151            state[0] + self.dt * k2[0],
152            state[1] + self.dt * k2[1],
153            state[2] + self.dt * k2[2],
154            state[3] + self.dt * k2[3],
155        ]
156    }
157
158    /// Advance one source-grid step atomically.
159    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
160        if !current.is_finite() {
161            return Err("current must be finite");
162        }
163        if !self.valid() {
164            return Err("NMDA state and parameters must satisfy the public bounds");
165        }
166        let held = self.refractory_remaining > 0.0;
167        let voltage = if held { self.v_reset } else { self.v };
168        let mut next = self.rk2_candidate(voltage, current);
169        let mut refractory = (self.refractory_remaining - self.dt).max(0.0);
170        let mut fired = 0;
171        if held {
172            next[0] = self.v_reset;
173        } else if next[0] >= self.v_threshold {
174            fired = 1;
175            next[0] = self.v_reset;
176            refractory = self.refractory_period;
177            next[1] += self.kinetic_scale * self.alpha_x;
178            next[3] += self.alpha_ca;
179        }
180        if !next.into_iter().chain([refractory]).all(f64::is_finite) {
181            return Err("NMDA candidate state became non-finite");
182        }
183        self.v = next[0].clamp(-120.0, 80.0);
184        self.x_nmda = next[1].max(0.0);
185        self.s_nmda = next[2].clamp(0.0, 1.0);
186        self.ca = next[3].max(0.0);
187        self.refractory_remaining = refractory;
188        Ok(fired)
189    }
190
191    /// Legacy fail-closed NetworkRunner surface.
192    pub fn step(&mut self, current: f64) -> i32 {
193        self.try_step(current).unwrap_or(0)
194    }
195
196    pub fn reset(&mut self) {
197        self.v = self.v_l;
198        self.x_nmda = 0.0;
199        self.s_nmda = 0.0;
200        self.ca = 0.0;
201        self.refractory_remaining = 0.0;
202    }
203}
204
205impl Default for NMDANeuron {
206    fn default() -> Self {
207        Self::new()
208    }
209}
210
211#[cfg(test)]
212mod tests {
213    use super::*;
214
215    #[test]
216    fn source_anchor_and_atomic_failure() {
217        let mut state = NMDANeuron::new();
218        assert_eq!(state.try_step(0.3), Ok(0));
219        assert!((state.v - -69.970_037_5).abs() < 1.0e-12);
220        let before = state.clone();
221        assert!(state.try_step(f64::NAN).is_err());
222        assert_eq!(state.v, before.v);
223        assert_eq!(state.s_nmda, before.s_nmda);
224    }
225}
sc_neurocore_engine/neurons/channels/nmda.rs

sc_neurocore_engine/neurons/channels/
nmda.rs
