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 — NMDA receptor-gated channel neuron
8
9use crate::neurons::biophysical::safe_rate;
10
11/// NMDA receptor neuron — WB base + NMDA-type glutamate receptor current.
12///
13/// NMDA receptors require both glutamate binding (modelled as input current)
14/// AND membrane depolarisation (Mg2+ block removal) for activation. The
15/// Mg2+ block is voltage-dependent: at rest (-65 mV) channels are blocked,
16/// but depolarisation to -40 mV relieves ~80% of the block.
17///
18/// Key mechanism for:
19/// - Coincidence detection: requires presynaptic (glutamate) + postsynaptic
20///   (depolarisation) signals simultaneously
21/// - Synaptic plasticity: Ca2+ influx through NMDA triggers LTP/LTD
22/// - Working memory: NMDA-mediated recurrent excitation sustains persistent
23///   activity in prefrontal cortex
24/// - Slow synaptic integration: rise ~10 ms, decay ~100 ms
25///
26/// Jahr & Stevens, J Neurosci 10:1830, 1990; Wang, Neuron 22:409, 1999.
27#[derive(Clone, Debug)]
28pub struct NMDANeuron {
29    pub v: f64,
30    pub h: f64,
31    pub n: f64,
32    pub s_nmda: f64, // NMDA synaptic variable (slow rise/decay)
33    pub g_na: f64,
34    pub g_k: f64,
35    pub g_nmda: f64, // NMDA conductance
36    pub g_l: f64,
37    pub e_na: f64,
38    pub e_k: f64,
39    pub e_nmda: f64, // NMDA reversal (0 mV, mixed cation)
40    pub e_l: f64,
41    pub c_m: f64,
42    pub phi: f64,
43    pub mg_conc: f64,   // Extracellular Mg2+ (mM), typically 1.0
44    pub tau_rise: f64,  // NMDA rise time (ms)
45    pub tau_decay: f64, // NMDA decay time (ms)
46    pub dt: f64,
47    pub v_threshold: f64,
48    pub gain: f64,
49}
50
51impl Default for NMDANeuron {
52    fn default() -> Self {
53        Self::new()
54    }
55}
56
57impl NMDANeuron {
58    pub fn new() -> Self {
59        Self {
60            v: -65.0,
61            h: 0.6,
62            n: 0.32,
63            s_nmda: 0.0,
64            g_na: 35.0,
65            g_k: 9.0,
66            g_nmda: 0.5,
67            g_l: 0.1,
68            e_na: 55.0,
69            e_k: -90.0,
70            e_nmda: 0.0, // Mixed cation reversal
71            e_l: -65.0,
72            c_m: 1.0,
73            phi: 5.0,
74            mg_conc: 1.0,
75            tau_rise: 10.0,
76            tau_decay: 100.0,
77            dt: 0.5,
78            v_threshold: -20.0,
79            gain: 1.0,
80        }
81    }
82
83    pub fn step(&mut self, current: f64) -> i32 {
84        let input = self.gain * current;
85        let sub_steps = 50;
86        let sub_dt = self.dt / sub_steps as f64;
87        let mut fired = 0i32;
88
89        // NMDA synaptic variable: driven by input (as proxy for glutamate)
90        let drive = if input > 0.0 {
91            input / (input + 5.0)
92        } else {
93            0.0
94        };
95        let ds = (drive - self.s_nmda)
96            / if drive > self.s_nmda {
97                self.tau_rise
98            } else {
99                self.tau_decay
100            };
101        self.s_nmda += self.dt * ds;
102        self.s_nmda = self.s_nmda.clamp(0.0, 1.0);
103
104        for _ in 0..sub_steps {
105            let v = self.v;
106
107            let alpha_m = safe_rate(0.1, 35.0, v, 10.0, 1.0);
108            let beta_m = 4.0 * (-(v + 60.0) / 18.0).exp();
109            let m_inf = alpha_m / (alpha_m + beta_m);
110
111            let alpha_h = 0.07 * (-(v + 58.0) / 20.0).exp();
112            let beta_h = 1.0 / (1.0 + (-(v + 28.0) / 10.0).exp());
113
114            let alpha_n = safe_rate(0.01, 34.0, v, 10.0, 0.1);
115            let beta_n = 0.125 * (-(v + 44.0) / 80.0).exp();
116
117            // Mg2+ block: B(V) = 1 / (1 + [Mg2+]/3.57 * exp(-0.062 * V))
118            // Jahr & Stevens 1990
119            let mg_block = 1.0 / (1.0 + (self.mg_conc / 3.57) * (-0.062 * v).exp());
120
121            self.h += sub_dt * self.phi * (alpha_h * (1.0 - self.h) - beta_h * self.h);
122            self.n += sub_dt * self.phi * (alpha_n * (1.0 - self.n) - beta_n * self.n);
123
124            let i_na = self.g_na * m_inf.powi(3) * self.h * (v - self.e_na);
125            let i_k = self.g_k * self.n.powi(4) * (v - self.e_k);
126            let i_nmda = self.g_nmda * self.s_nmda * mg_block * (v - self.e_nmda);
127            let i_l = self.g_l * (v - self.e_l);
128
129            let dv = (-i_na - i_k - i_nmda - i_l + input) / self.c_m;
130            self.v += sub_dt * dv;
131
132            if self.v >= self.v_threshold {
133                fired = 1;
134                self.v = -65.0;
135            }
136        }
137
138        self.v = self.v.clamp(-100.0, 60.0);
139        if !self.v.is_finite() {
140            self.v = -65.0;
141            self.h = 0.6;
142            self.n = 0.32;
143        }
144        if !self.s_nmda.is_finite() {
145            self.s_nmda = 0.0;
146        }
147        self.h = self.h.clamp(0.0, 1.0);
148        self.n = self.n.clamp(0.0, 1.0);
149
150        fired
151    }
152
153    pub fn reset(&mut self) {
154        *self = Self::new();
155    }
156}
157
158#[cfg(test)]
159mod tests {
160    use super::*;
161
162    // -- NMDA Neuron tests --
163
164    #[test]
165    fn nmda_fires_with_input() {
166        let mut n = NMDANeuron::new();
167        let mut spikes = 0;
168        for _ in 0..2_000 {
169            spikes += n.step(3.0);
170        }
171        assert!(spikes > 5, "NMDA neuron must fire with input, got {spikes}");
172    }
173
174    #[test]
175    fn nmda_silent_without_input() {
176        let mut n = NMDANeuron::new();
177        let mut spikes = 0;
178        for _ in 0..10_000 {
179            spikes += n.step(0.0);
180        }
181        assert_eq!(
182            spikes, 0,
183            "NMDA neuron must be silent without input, got {spikes}"
184        );
185    }
186
187    #[test]
188    fn nmda_mg_block_at_rest() {
189        // At -65 mV: B = 1/(1 + 1/3.57 * exp(0.062*65)) = 1/(1 + 0.28 * 56.3) = 1/16.8 = 0.06
190        let n = NMDANeuron::new();
191        let mg_block = 1.0 / (1.0 + (n.mg_conc / 3.57) * (-0.062 * n.v).exp());
192        assert!(
193            mg_block < 0.1,
194            "Mg2+ block must be strong at rest, B={mg_block}"
195        );
196    }
197
198    #[test]
199    fn nmda_mg_relief_at_depolarised() {
200        // At -20 mV: B = 1/(1 + 0.28 * exp(0.062*20)) = 1/(1 + 0.28*3.45) = 1/1.97 = 0.51
201        let mg_block = 1.0 / (1.0 + (1.0 / 3.57) * (-0.062 * (-20.0_f64)).exp());
202        assert!(
203            mg_block > 0.4,
204            "Mg2+ block must be relieved at -20 mV, B={mg_block}"
205        );
206    }
207
208    #[test]
209    fn nmda_s_builds_with_input() {
210        let mut n = NMDANeuron::new();
211        assert_eq!(n.s_nmda, 0.0);
212        for _ in 0..2000 {
213            n.step(5.0);
214        }
215        assert!(
216            n.s_nmda > 0.0,
217            "s_nmda must build with input, s={}",
218            n.s_nmda
219        );
220    }
221
222    #[test]
223    fn nmda_s_decays_without_input() {
224        let mut n = NMDANeuron::new();
225        // Build up s
226        for _ in 0..2000 {
227            n.step(5.0);
228        }
229        let s_peak = n.s_nmda;
230        // Remove input
231        for _ in 0..2000 {
232            n.step(0.0);
233        }
234        assert!(n.s_nmda < s_peak, "s_nmda must decay after input removal");
235    }
236
237    #[test]
238    fn nmda_zero_mg_increases_current() {
239        // Without Mg2+ block, NMDA should contribute more
240        let mut with_mg = NMDANeuron::new();
241        let mut no_mg = NMDANeuron::new();
242        no_mg.mg_conc = 0.0;
243
244        let input = 2.0;
245        let mut spikes_mg = 0;
246        let mut spikes_no = 0;
247        for _ in 0..10_000 {
248            spikes_mg += with_mg.step(input);
249            spikes_no += no_mg.step(input);
250        }
251        assert!(
252            spikes_no >= spikes_mg,
253            "No Mg2+ should increase NMDA current: no_mg={spikes_no} vs mg={spikes_mg}"
254        );
255    }
256
257    #[test]
258    fn nmda_negative_input_no_crash() {
259        let mut n = NMDANeuron::new();
260        for _ in 0..10_000 {
261            n.step(-100.0);
262        }
263        assert!(n.v.is_finite());
264    }
265
266    #[test]
267    fn nmda_nan_input_stays_finite() {
268        let mut n = NMDANeuron::new();
269        n.step(f64::NAN);
270        assert!(n.v.is_finite());
271    }
272
273    #[test]
274    fn nmda_extreme_input_bounded() {
275        let mut n = NMDANeuron::new();
276        for _ in 0..1000 {
277            n.step(1e6);
278        }
279        assert!(n.v.is_finite() && n.v <= 60.0);
280    }
281
282    #[test]
283    fn nmda_reset_clears_state() {
284        let mut n = NMDANeuron::new();
285        for _ in 0..1000 {
286            n.step(10.0);
287        }
288        n.reset();
289        assert_eq!(n.v, -65.0);
290        assert_eq!(n.s_nmda, 0.0);
291    }
292
293    #[test]
294    fn nmda_performance_1k_steps() {
295        let start = std::time::Instant::now();
296        let mut n = NMDANeuron::new();
297        for _ in 0..1_000 {
298            std::hint::black_box(n.step(3.0));
299        }
300        let elapsed = start.elapsed();
301        assert!(
302            elapsed.as_millis() < 200,
303            "1k steps must complete in <200ms"
304        );
305    }
306}
sc_neurocore_engine/neurons/channels/nmda.rs

sc_neurocore_engine/neurons/channels/
nmda.rs
