sc_neurocore_engine/neurons/rate/
compte_wm.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 — Compte et al. 2000 pyramidal cell
8
9//! Source-bounded Compte pyramidal-cell and incoming channel dynamics.
10
11const V_MIN: f64 = -200.0;
12const V_MAX: f64 = 100.0;
13const GATE_MAX: f64 = 1.0e6;
14
15/// Complete Compte pyramidal-cell state and source control-set configuration.
16///
17/// Conductances are microSiemens, voltages millivolts, current nanoamps,
18/// capacitance nanofarads, and time milliseconds. Incoming event flags map to
19/// recurrent NMDA, external AMPA, and inhibitory GABAA pathways respectively.
20#[derive(Clone, Debug)]
21pub struct CompteWMNeuron {
22    /// Membrane potential in mV.
23    pub v: f64,
24    /// External AMPA open-channel accumulator.
25    pub s_ampa: f64,
26    /// Recurrent NMDA open fraction.
27    pub s_nmda: f64,
28    /// Recurrent NMDA rise precursor.
29    pub x_nmda: f64,
30    /// Incoming GABAA open-channel accumulator.
31    pub s_gaba: f64,
32    /// Remaining absolute refractory duration in ms.
33    pub ref_remaining: f64,
34    /// Leak conductance in microSiemens.
35    pub g_l: f64,
36    /// External pyramidal AMPA conductance in microSiemens.
37    pub g_ampa: f64,
38    /// Recurrent pyramidal NMDA conductance in microSiemens.
39    pub g_nmda: f64,
40    /// Interneuron-to-pyramidal GABAA conductance in microSiemens.
41    pub g_gaba: f64,
42    /// Leak reversal in mV.
43    pub e_l: f64,
44    /// Excitatory reversal in mV.
45    pub e_exc: f64,
46    /// Inhibitory reversal in mV.
47    pub e_inh: f64,
48    /// Membrane capacitance in nF.
49    pub c_m: f64,
50    /// Extracellular magnesium concentration in mM.
51    pub mg: f64,
52    /// AMPA decay constant in ms.
53    pub tau_ampa: f64,
54    /// NMDA decay constant in ms.
55    pub tau_nmda: f64,
56    /// NMDA rise-precursor decay constant in ms.
57    pub tau_x: f64,
58    /// GABAA decay constant in ms.
59    pub tau_gaba: f64,
60    /// NMDA saturation rate in inverse ms.
61    pub alpha_nmda: f64,
62    /// Sampled firing threshold in mV.
63    pub v_threshold: f64,
64    /// Post-spike and refractory voltage in mV.
65    pub v_reset: f64,
66    /// Absolute refractory duration in ms.
67    pub tau_ref: f64,
68    /// Midpoint-RK2 integration step in ms.
69    pub dt: f64,
70}
71
72impl CompteWMNeuron {
73    /// Construct the Compte (2000) source control-set pyramidal defaults.
74    #[must_use]
75    pub fn new() -> Self {
76        Self {
77            v: -70.0,
78            s_ampa: 0.0,
79            s_nmda: 0.0,
80            x_nmda: 0.0,
81            s_gaba: 0.0,
82            ref_remaining: 0.0,
83            g_l: 0.025,
84            g_ampa: 0.0031,
85            g_nmda: 0.000_381,
86            g_gaba: 0.001_336,
87            e_l: -70.0,
88            e_exc: 0.0,
89            e_inh: -70.0,
90            c_m: 0.5,
91            mg: 1.0,
92            tau_ampa: 2.0,
93            tau_nmda: 100.0,
94            tau_x: 2.0,
95            tau_gaba: 10.0,
96            alpha_nmda: 0.5,
97            v_threshold: -50.0,
98            v_reset: -60.0,
99            tau_ref: 2.0,
100            dt: 0.02,
101        }
102    }
103
104    /// Return whether all mutable state and configuration invariants hold.
105    #[must_use]
106    pub fn validate(&self) -> bool {
107        let finite = [
108            self.v,
109            self.s_ampa,
110            self.s_nmda,
111            self.x_nmda,
112            self.s_gaba,
113            self.ref_remaining,
114            self.g_l,
115            self.g_ampa,
116            self.g_nmda,
117            self.g_gaba,
118            self.e_l,
119            self.e_exc,
120            self.e_inh,
121            self.c_m,
122            self.mg,
123            self.tau_ampa,
124            self.tau_nmda,
125            self.tau_x,
126            self.tau_gaba,
127            self.alpha_nmda,
128            self.v_threshold,
129            self.v_reset,
130            self.tau_ref,
131            self.dt,
132        ]
133        .iter()
134        .all(|value| value.is_finite());
135        finite
136            && (V_MIN..=V_MAX).contains(&self.v)
137            && (V_MIN..=V_MAX).contains(&self.v_reset)
138            && [self.s_ampa, self.x_nmda, self.s_gaba]
139                .iter()
140                .all(|value| (0.0..=GATE_MAX).contains(value))
141            && (0.0..=1.0).contains(&self.s_nmda)
142            && self.ref_remaining >= 0.0
143            && [
144                self.g_l,
145                self.g_ampa,
146                self.g_nmda,
147                self.g_gaba,
148                self.mg,
149                self.alpha_nmda,
150            ]
151            .iter()
152            .all(|value| *value >= 0.0)
153            && [
154                self.c_m,
155                self.tau_ampa,
156                self.tau_nmda,
157                self.tau_x,
158                self.tau_gaba,
159                self.tau_ref,
160                self.dt,
161            ]
162            .iter()
163            .all(|value| *value > 0.0)
164    }
165
166    fn mg_block(&self, v: f64) -> Option<f64> {
167        let exponent = -0.062 * v;
168        let block = if exponent > 700.0 {
169            0.0
170        } else {
171            1.0 / (1.0 + self.mg / 3.57 * exponent.exp())
172        };
173        (block.is_finite() && (0.0..=1.0).contains(&block)).then_some(block)
174    }
175
176    fn derivatives(
177        &self,
178        state: [f64; 5],
179        current: f64,
180        membrane_active: bool,
181    ) -> Option<[f64; 5]> {
182        let [v, s_ampa, s_nmda, x_nmda, s_gaba] = state;
183        let d_ampa = -s_ampa / self.tau_ampa;
184        let d_nmda = -s_nmda / self.tau_nmda + self.alpha_nmda * x_nmda * (1.0 - s_nmda);
185        let d_x = -x_nmda / self.tau_x;
186        let d_gaba = -s_gaba / self.tau_gaba;
187        let d_v = if membrane_active {
188            let i_l = self.g_l * (v - self.e_l);
189            let i_ampa = self.g_ampa * s_ampa * (v - self.e_exc);
190            let i_nmda = self.g_nmda * self.mg_block(v)? * s_nmda * (v - self.e_exc);
191            let i_gaba = self.g_gaba * s_gaba * (v - self.e_inh);
192            (-i_l - i_ampa - i_nmda - i_gaba + current) / self.c_m
193        } else {
194            0.0
195        };
196        let result = [d_v, d_ampa, d_nmda, d_x, d_gaba];
197        result
198            .iter()
199            .all(|value| value.is_finite())
200            .then_some(result)
201    }
202
203    /// Advance one atomic source-level midpoint-RK2 step.
204    ///
205    /// The three flags apply recurrent-NMDA, external-AMPA, and
206    /// inhibitory-GABAA event jumps before the continuous flow. Threshold
207    /// detection is sampled; within-step firing-time interpolation is not
208    /// claimed. An error leaves all public state unchanged.
209    pub fn step_events(
210        &mut self,
211        current: f64,
212        recurrent_event: bool,
213        external_event: bool,
214        inhibitory_event: bool,
215    ) -> Result<i32, &'static str> {
216        if !self.validate() || !current.is_finite() {
217            return Err("invalid Compte state, configuration, or current");
218        }
219        let initial = [
220            self.v,
221            self.s_ampa + if external_event { 1.0 } else { 0.0 },
222            self.s_nmda,
223            self.x_nmda + if recurrent_event { 1.0 } else { 0.0 },
224            self.s_gaba + if inhibitory_event { 1.0 } else { 0.0 },
225        ];
226        if !initial[1..]
227            .iter()
228            .all(|value| value.is_finite() && (0.0..=GATE_MAX).contains(value))
229        {
230            return Err("Compte event candidate outside gate envelope");
231        }
232        let active = self.ref_remaining <= 0.0;
233        let k1 = self
234            .derivatives(initial, current, active)
235            .ok_or("non-finite Compte RK2 first stage")?;
236        let midpoint = std::array::from_fn(|index| initial[index] + 0.5 * self.dt * k1[index]);
237        let k2 = self
238            .derivatives(midpoint, current, active)
239            .ok_or("non-finite Compte RK2 midpoint stage")?;
240        let mut candidate: [f64; 5] =
241            std::array::from_fn(|index| initial[index] + self.dt * k2[index]);
242        if !candidate.iter().all(|value| value.is_finite())
243            || !(V_MIN..=V_MAX).contains(&candidate[0])
244            || !candidate[1..]
245                .iter()
246                .all(|value| (0.0..=GATE_MAX).contains(value))
247            || candidate[2] > 1.0
248        {
249            return Err("Compte RK2 candidate outside safety envelope");
250        }
251        let mut event = 0;
252        let mut ref_remaining = (self.ref_remaining - self.dt).max(0.0);
253        if !active {
254            candidate[0] = self.v_reset;
255        } else if candidate[0] >= self.v_threshold {
256            candidate[0] = self.v_reset;
257            ref_remaining = self.tau_ref;
258            event = 1;
259        }
260        self.v = candidate[0];
261        self.s_ampa = candidate[1];
262        self.s_nmda = candidate[2];
263        self.x_nmda = candidate[3];
264        self.s_gaba = candidate[4];
265        self.ref_remaining = ref_remaining;
266        Ok(event)
267    }
268
269    /// Retain the generic two-argument engine adapter.
270    ///
271    /// The boolean retains its historical meaning as a recurrent excitatory
272    /// NMDA event. Invalid state returns zero without mutation; callers that
273    /// require explicit errors use the complete step_events method.
274    pub fn step(&mut self, current: f64, recurrent_event: bool) -> i32 {
275        self.step_events(current, recurrent_event, false, false)
276            .unwrap_or(0)
277    }
278
279    /// Reset all dynamic state while preserving configuration.
280    pub fn reset(&mut self) {
281        self.v = self.e_l;
282        self.s_ampa = 0.0;
283        self.s_nmda = 0.0;
284        self.x_nmda = 0.0;
285        self.s_gaba = 0.0;
286        self.ref_remaining = 0.0;
287    }
288
289    /// Return the complete dynamic state in public trace order.
290    #[must_use]
291    pub fn get_state(&self) -> [f64; 6] {
292        [
293            self.v,
294            self.s_ampa,
295            self.s_nmda,
296            self.x_nmda,
297            self.s_gaba,
298            self.ref_remaining,
299        ]
300    }
301}
302
303impl Default for CompteWMNeuron {
304    fn default() -> Self {
305        Self::new()
306    }
307}
308
309#[cfg(test)]
310mod tests {
311    use super::*;
312
313    #[test]
314    fn event_pathways_are_separate() {
315        let mut n = CompteWMNeuron::new();
316        assert_eq!(n.step_events(0.0, true, false, false), Ok(0));
317        assert_eq!(n.s_ampa, 0.0);
318        assert!(n.s_nmda > 0.0 && n.x_nmda > 0.0);
319        assert_eq!(n.s_gaba, 0.0);
320    }
321
322    #[test]
323    fn failure_is_atomic() {
324        let mut n = CompteWMNeuron::new();
325        let before = n.get_state();
326        assert!(n.step_events(f64::NAN, false, false, false).is_err());
327        assert_eq!(n.get_state(), before);
328    }
329
330    #[test]
331    fn reset_preserves_configuration() {
332        let mut n = CompteWMNeuron::new();
333        n.dt = 0.01;
334        n.step_events(1.0, true, true, true).unwrap();
335        n.reset();
336        assert_eq!(n.get_state(), [-70.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
337        assert_eq!(n.dt, 0.01);
338    }
339}
sc_neurocore_engine/neurons/rate/compte_wm.rs

sc_neurocore_engine/neurons/rate/
compte_wm.rs
