sc_neurocore_engine/neurons/
jansen_rit.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 — Jansen and Rit 1995 cortical-column neural mass
8
9//! Equation-(6) dynamics and an atomic explicit-Euler batch implementation.
10
11/// One Jansen–Rit cortical column with states `[y0, y1, y2, y3, y4, y5]`.
12#[derive(Clone, Debug)]
13pub struct JansenRitUnit {
14    pub y: [f64; 6],
15    pub a_exc: f64,
16    pub b_exc: f64,
17    pub a_rate: f64,
18    pub b_rate: f64,
19    pub c: f64,
20    pub e0: f64,
21    pub v0: f64,
22    pub r: f64,
23    pub dt: f64,
24}
25
26impl JansenRitUnit {
27    /// Construct the published parameter set with a 0.1 ms Euler step.
28    pub fn new() -> Self {
29        Self {
30            y: [0.0; 6],
31            a_exc: 3.25,
32            b_exc: 22.0,
33            a_rate: 100.0,
34            b_rate: 50.0,
35            c: 135.0,
36            e0: 2.5,
37            v0: 6.0,
38            r: 0.56,
39            dt: 0.0001,
40        }
41    }
42
43    /// Construct and validate one configured state.
44    #[allow(clippy::too_many_arguments)]
45    pub fn with_parameters(
46        y0: f64,
47        y3: f64,
48        y1: f64,
49        y4: f64,
50        y2: f64,
51        y5: f64,
52        a_exc: f64,
53        b_exc: f64,
54        a_rate: f64,
55        b_rate: f64,
56        c: f64,
57        e0: f64,
58        v0: f64,
59        r: f64,
60        dt: f64,
61    ) -> Result<Self, String> {
62        let unit = Self {
63            y: [y0, y1, y2, y3, y4, y5],
64            a_exc,
65            b_exc,
66            a_rate,
67            b_rate,
68            c,
69            e0,
70            v0,
71            r,
72            dt,
73        };
74        unit.validate()?;
75        Ok(unit)
76    }
77
78    fn validate(&self) -> Result<(), String> {
79        let parameters = [
80            self.a_exc,
81            self.b_exc,
82            self.a_rate,
83            self.b_rate,
84            self.c,
85            self.e0,
86            self.v0,
87            self.r,
88            self.dt,
89        ];
90        if !self
91            .y
92            .iter()
93            .chain(parameters.iter())
94            .all(|value| value.is_finite())
95        {
96            return Err("Jansen–Rit state and parameters must be finite".into());
97        }
98        if self.a_exc <= 0.0
99            || self.b_exc <= 0.0
100            || self.a_rate <= 0.0
101            || self.b_rate <= 0.0
102            || self.e0 <= 0.0
103            || self.r <= 0.0
104            || self.dt <= 0.0
105        {
106            return Err(
107                "Jansen–Rit gains, rates, sigmoid scale, slope, and dt must be positive".into(),
108            );
109        }
110        if self.c < 0.0 {
111            return Err("Jansen–Rit connectivity must be non-negative".into());
112        }
113        Ok(())
114    }
115
116    #[inline]
117    fn sigmoid(&self, voltage: f64) -> Result<f64, String> {
118        if !voltage.is_finite() {
119            return Err("Jansen–Rit sigmoid input must be finite".into());
120        }
121        let exponent = self.r * (self.v0 - voltage);
122        let response = if exponent >= 0.0 {
123            let exp_neg = (-exponent).exp();
124            2.0 * self.e0 * exp_neg / (1.0 + exp_neg)
125        } else {
126            2.0 * self.e0 / (1.0 + exponent.exp())
127        };
128        if !response.is_finite() {
129            return Err("Jansen–Rit sigmoid response must be finite".into());
130        }
131        Ok(response)
132    }
133
134    /// Advance one equation-(6) Euler step and return post-update `y1 - y2`.
135    pub fn step(&mut self, p_ext: f64) -> Result<f64, String> {
136        self.validate()?;
137        if !p_ext.is_finite() {
138            return Err("Jansen–Rit external drive must be finite".into());
139        }
140        let c1 = self.c;
141        let c2 = 0.8 * c1;
142        let c3 = 0.25 * c1;
143        let c4 = 0.25 * c1;
144        let s_pyramidal = self.sigmoid(self.y[1] - self.y[2])?;
145        let s_excitatory = self.sigmoid(c1 * self.y[0])?;
146        let s_inhibitory = self.sigmoid(c3 * self.y[0])?;
147        let derivatives = [
148            self.y[3],
149            self.y[4],
150            self.y[5],
151            self.a_exc * self.a_rate * s_pyramidal
152                - 2.0 * self.a_rate * self.y[3]
153                - self.a_rate.powi(2) * self.y[0],
154            self.a_exc * self.a_rate * (p_ext + c2 * s_excitatory)
155                - 2.0 * self.a_rate * self.y[4]
156                - self.a_rate.powi(2) * self.y[1],
157            self.b_exc * self.b_rate * c4 * s_inhibitory
158                - 2.0 * self.b_rate * self.y[5]
159                - self.b_rate.powi(2) * self.y[2],
160        ];
161        let mut candidate = self.y;
162        for (next, derivative) in candidate.iter_mut().zip(derivatives) {
163            *next += self.dt * derivative;
164        }
165        if !candidate.iter().all(|value| value.is_finite()) {
166            return Err("Jansen–Rit candidate state must remain finite".into());
167        }
168        self.y = candidate;
169        Ok(self.y[1] - self.y[2])
170    }
171
172    /// Restore all dynamic states while preserving parameters.
173    pub fn reset(&mut self) {
174        self.y = [0.0; 6];
175    }
176}
177
178impl Default for JansenRitUnit {
179    fn default() -> Self {
180        Self::new()
181    }
182}
183
184/// Per-step state and EEG traces returned by the batch implementation.
185pub struct JansenRitTrace {
186    pub y0: Vec<f64>,
187    pub y3: Vec<f64>,
188    pub y1: Vec<f64>,
189    pub y4: Vec<f64>,
190    pub y2: Vec<f64>,
191    pub y5: Vec<f64>,
192    pub eeg: Vec<f64>,
193    pub final_state: [f64; 6],
194}
195
196/// Simulate a complete external-drive batch.
197#[allow(clippy::too_many_arguments)]
198pub fn simulate(
199    y0: f64,
200    y3: f64,
201    y1: f64,
202    y4: f64,
203    y2: f64,
204    y5: f64,
205    a_exc: f64,
206    b_exc: f64,
207    a_rate: f64,
208    b_rate: f64,
209    c: f64,
210    e0: f64,
211    v0: f64,
212    r: f64,
213    dt: f64,
214    p_ext: &[f64],
215) -> Result<JansenRitTrace, String> {
216    let mut unit = JansenRitUnit::with_parameters(
217        y0, y3, y1, y4, y2, y5, a_exc, b_exc, a_rate, b_rate, c, e0, v0, r, dt,
218    )?;
219    let mut trace = JansenRitTrace {
220        y0: Vec::with_capacity(p_ext.len()),
221        y3: Vec::with_capacity(p_ext.len()),
222        y1: Vec::with_capacity(p_ext.len()),
223        y4: Vec::with_capacity(p_ext.len()),
224        y2: Vec::with_capacity(p_ext.len()),
225        y5: Vec::with_capacity(p_ext.len()),
226        eeg: Vec::with_capacity(p_ext.len()),
227        final_state: unit.y,
228    };
229    for drive in p_ext {
230        let eeg = unit.step(*drive)?;
231        trace.y0.push(unit.y[0]);
232        trace.y3.push(unit.y[3]);
233        trace.y1.push(unit.y[1]);
234        trace.y4.push(unit.y[4]);
235        trace.y2.push(unit.y[2]);
236        trace.y5.push(unit.y[5]);
237        trace.eeg.push(eeg);
238    }
239    trace.final_state = unit.y;
240    Ok(trace)
241}
242
243#[cfg(test)]
244mod tests {
245    use super::*;
246
247    #[test]
248    fn one_step_uses_c1_inside_and_c2_outside_excitatory_sigmoid() {
249        let mut unit = JansenRitUnit::with_parameters(
250            0.1, 0.2, 0.3, -0.4, -0.1, 0.5, 3.25, 22.0, 100.0, 50.0, 135.0, 2.5, 6.0, 0.56, 0.0001,
251        )
252        .unwrap();
253        let old = unit.y;
254        let se = unit.sigmoid(135.0 * old[0]).unwrap();
255        let expected_y4 = old[4]
256            + 0.0001
257                * (3.25 * 100.0 * (220.0 + 0.8 * 135.0 * se)
258                    - 2.0 * 100.0 * old[4]
259                    - 100.0_f64.powi(2) * old[1]);
260        unit.step(220.0).unwrap();
261        assert_eq!(unit.y[4], expected_y4);
262    }
263
264    #[test]
265    fn batch_matches_scalar_and_preserves_empty_initial_state() {
266        let empty = simulate(
267            0.1,
268            0.2,
269            0.3,
270            -0.4,
271            -0.1,
272            0.5,
273            3.25,
274            22.0,
275            100.0,
276            50.0,
277            135.0,
278            2.5,
279            6.0,
280            0.56,
281            0.0001,
282            &[],
283        )
284        .unwrap();
285        assert!(empty.eeg.is_empty());
286        assert_eq!(empty.final_state, [0.1, 0.3, -0.1, 0.2, -0.4, 0.5]);
287
288        let drives = [120.0, 220.0, 320.0];
289        let batch = simulate(
290            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.25, 22.0, 100.0, 50.0, 135.0, 2.5, 6.0, 0.56, 0.0001,
291            &drives,
292        )
293        .unwrap();
294        let mut scalar = JansenRitUnit::new();
295        for drive in drives {
296            scalar.step(drive).unwrap();
297        }
298        assert_eq!(batch.final_state, scalar.y);
299    }
300
301    #[test]
302    fn invalid_input_does_not_mutate_state() {
303        let mut unit = JansenRitUnit::new();
304        let before = unit.y;
305        assert!(unit.step(f64::NAN).is_err());
306        assert_eq!(unit.y, before);
307    }
308}
sc_neurocore_engine/neurons/jansen_rit.rs

sc_neurocore_engine/neurons/
jansen_rit.rs
