sc_neurocore_engine/
wong_wang.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 — Wong and Wang 2006 reduced decision circuit
8
9//! Explicit-Euler NMDA dynamics with the published AMPA Ornstein-Uhlenbeck
10//! current noise.  The deterministic sample-taking boundary is shared by the
11//! scalar engine class and the batch accelerator.
12
13use rand::RngExt;
14use rand::SeedableRng;
15use rand_xoshiro::Xoshiro256PlusPlus;
16
17const A: f64 = 270.0;
18const B: f64 = 108.0;
19const D: f64 = 0.154;
20
21/// Four dynamic states of one reduced Wong-Wang decision circuit.
22#[derive(Clone, Debug)]
23pub struct WongWangUnit {
24    pub s1: f64,
25    pub s2: f64,
26    pub noise1: f64,
27    pub noise2: f64,
28    pub tau_s: f64,
29    pub tau_ampa: f64,
30    pub gamma: f64,
31    pub j_n: f64,
32    pub j_cross: f64,
33    pub i_0: f64,
34    pub sigma: f64,
35    pub dt: f64,
36    rng: Xoshiro256PlusPlus,
37}
38
39impl WongWangUnit {
40    /// Construct the published parameter set with a reproducible RNG stream.
41    pub fn new(seed: u64) -> Self {
42        Self {
43            s1: 0.1,
44            s2: 0.1,
45            noise1: 0.0,
46            noise2: 0.0,
47            tau_s: 0.1,
48            tau_ampa: 0.002,
49            gamma: 0.641,
50            j_n: 0.2609,
51            j_cross: 0.0497,
52            i_0: 0.3255,
53            sigma: 0.02,
54            dt: 0.0001,
55            rng: Xoshiro256PlusPlus::seed_from_u64(seed),
56        }
57    }
58
59    /// Construct one configured state while retaining the supplied RNG seed.
60    #[allow(clippy::too_many_arguments)]
61    pub fn with_parameters(
62        s1: f64,
63        s2: f64,
64        noise1: f64,
65        noise2: f64,
66        tau_s: f64,
67        tau_ampa: f64,
68        gamma: f64,
69        j_n: f64,
70        j_cross: f64,
71        i_0: f64,
72        sigma: f64,
73        dt: f64,
74        seed: u64,
75    ) -> Result<Self, String> {
76        let unit = Self {
77            s1,
78            s2,
79            noise1,
80            noise2,
81            tau_s,
82            tau_ampa,
83            gamma,
84            j_n,
85            j_cross,
86            i_0,
87            sigma,
88            dt,
89            rng: Xoshiro256PlusPlus::seed_from_u64(seed),
90        };
91        unit.validate()?;
92        Ok(unit)
93    }
94
95    fn validate(&self) -> Result<(), String> {
96        let finite = [
97            self.s1,
98            self.s2,
99            self.noise1,
100            self.noise2,
101            self.tau_s,
102            self.tau_ampa,
103            self.gamma,
104            self.j_n,
105            self.j_cross,
106            self.i_0,
107            self.sigma,
108            self.dt,
109        ];
110        if !finite.iter().all(|value| value.is_finite()) {
111            return Err("Wong-Wang state and parameters must be finite".into());
112        }
113        if !(0.0..=1.0).contains(&self.s1) || !(0.0..=1.0).contains(&self.s2) {
114            return Err("Wong-Wang gating state must remain in [0, 1]".into());
115        }
116        if self.tau_s <= 0.0 || self.tau_ampa <= 0.0 || self.gamma <= 0.0 || self.dt <= 0.0 {
117            return Err("Wong-Wang time constants, gamma, and dt must be positive".into());
118        }
119        if self.j_n < 0.0 || self.j_cross < 0.0 || self.sigma < 0.0 {
120            return Err("Wong-Wang couplings and sigma must be non-negative".into());
121        }
122        Ok(())
123    }
124
125    #[inline]
126    fn phi(i_syn: f64) -> Result<f64, String> {
127        if !i_syn.is_finite() {
128            return Err("Wong-Wang synaptic current must be finite".into());
129        }
130        let x = A * i_syn - B;
131        let scaled = -D * x;
132        let response = if scaled > 700.0 {
133            0.0
134        } else if x.abs() < 1.0e-7 {
135            1.0 / D
136        } else {
137            x / -scaled.exp_m1()
138        };
139        if !response.is_finite() {
140            return Err("Wong-Wang transfer response must be finite".into());
141        }
142        Ok(response.max(0.0))
143    }
144
145    /// Advance one Euler/OU update from externally supplied normal samples.
146    pub fn step_with_gaussian_samples(
147        &mut self,
148        stim1: f64,
149        stim2: f64,
150        xi1: f64,
151        xi2: f64,
152    ) -> Result<(f64, f64), String> {
153        self.validate()?;
154        if ![stim1, stim2, xi1, xi2]
155            .iter()
156            .all(|value| value.is_finite())
157        {
158            return Err("Wong-Wang stimuli and Gaussian samples must be finite".into());
159        }
160        let current1 = self.j_n * self.s1 - self.j_cross * self.s2 + self.i_0 + stim1 + self.noise1;
161        let current2 = self.j_n * self.s2 - self.j_cross * self.s1 + self.i_0 + stim2 + self.noise2;
162        let rate1 = Self::phi(current1)?;
163        let rate2 = Self::phi(current2)?;
164        let ds1 = -self.s1 / self.tau_s + (1.0 - self.s1) * self.gamma * rate1;
165        let ds2 = -self.s2 / self.tau_s + (1.0 - self.s2) * self.gamma * rate2;
166        let noise_scale = (self.dt / self.tau_ampa).sqrt() * self.sigma;
167        let next_s1 = self.s1 + self.dt * ds1;
168        let next_s2 = self.s2 + self.dt * ds2;
169        let next_noise1 = self.noise1 - (self.dt / self.tau_ampa) * self.noise1 + noise_scale * xi1;
170        let next_noise2 = self.noise2 - (self.dt / self.tau_ampa) * self.noise2 + noise_scale * xi2;
171        if ![next_s1, next_s2, next_noise1, next_noise2]
172            .iter()
173            .all(|value| value.is_finite())
174        {
175            return Err("Wong-Wang candidate state must remain finite".into());
176        }
177        if !(0.0..=1.0).contains(&next_s1) || !(0.0..=1.0).contains(&next_s2) {
178            return Err("Wong-Wang candidate gating state left [0, 1]".into());
179        }
180        self.s1 = next_s1;
181        self.s2 = next_s2;
182        self.noise1 = next_noise1;
183        self.noise2 = next_noise2;
184        Ok((rate1, rate2))
185    }
186
187    /// Advance one stochastic update using the internal seeded RNG.
188    pub fn step(&mut self, stim1: f64, stim2: f64) -> Result<(f64, f64), String> {
189        let xi1 = self.randn();
190        let xi2 = self.randn();
191        self.step_with_gaussian_samples(stim1, stim2, xi1, xi2)
192    }
193
194    fn randn(&mut self) -> f64 {
195        let u1 = self.rng.random::<f64>().max(f64::MIN_POSITIVE);
196        let u2 = self.rng.random::<f64>();
197        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
198    }
199
200    /// Restore only dynamic state variables.
201    pub fn reset(&mut self) {
202        self.s1 = 0.1;
203        self.s2 = 0.1;
204        self.noise1 = 0.0;
205        self.noise2 = 0.0;
206    }
207}
208
209/// Per-step state and rate traces returned by the batch implementation.
210pub struct WongWangTrace {
211    pub s1: Vec<f64>,
212    pub s2: Vec<f64>,
213    pub noise1: Vec<f64>,
214    pub noise2: Vec<f64>,
215    pub rate1: Vec<f64>,
216    pub rate2: Vec<f64>,
217    pub final_s1: f64,
218    pub final_s2: f64,
219    pub final_noise1: f64,
220    pub final_noise2: f64,
221}
222
223/// Simulate a complete deterministic-sample batch.
224#[allow(clippy::too_many_arguments)]
225pub fn simulate(
226    s1: f64,
227    s2: f64,
228    noise1: f64,
229    noise2: f64,
230    tau_s: f64,
231    tau_ampa: f64,
232    gamma: f64,
233    j_n: f64,
234    j_cross: f64,
235    i_0: f64,
236    sigma: f64,
237    dt: f64,
238    stim1: &[f64],
239    stim2: &[f64],
240    xi: &[f64],
241) -> Result<WongWangTrace, String> {
242    let n_steps = stim1.len();
243    if stim2.len() != n_steps {
244        return Err(format!(
245            "stim1 and stim2 length mismatch: {n_steps} vs {}",
246            stim2.len()
247        ));
248    }
249    if xi.len() != 2 * n_steps {
250        return Err(format!(
251            "xi length must be 2 * n_steps ({}): got {}",
252            2 * n_steps,
253            xi.len()
254        ));
255    }
256    let mut unit = WongWangUnit::with_parameters(
257        s1, s2, noise1, noise2, tau_s, tau_ampa, gamma, j_n, j_cross, i_0, sigma, dt, 0,
258    )?;
259    let mut trace = WongWangTrace {
260        s1: Vec::with_capacity(n_steps),
261        s2: Vec::with_capacity(n_steps),
262        noise1: Vec::with_capacity(n_steps),
263        noise2: Vec::with_capacity(n_steps),
264        rate1: Vec::with_capacity(n_steps),
265        rate2: Vec::with_capacity(n_steps),
266        final_s1: s1,
267        final_s2: s2,
268        final_noise1: noise1,
269        final_noise2: noise2,
270    };
271    for step in 0..n_steps {
272        let (rate1, rate2) = unit.step_with_gaussian_samples(
273            stim1[step],
274            stim2[step],
275            xi[2 * step],
276            xi[2 * step + 1],
277        )?;
278        trace.s1.push(unit.s1);
279        trace.s2.push(unit.s2);
280        trace.noise1.push(unit.noise1);
281        trace.noise2.push(unit.noise2);
282        trace.rate1.push(rate1);
283        trace.rate2.push(rate2);
284    }
285    trace.final_s1 = unit.s1;
286    trace.final_s2 = unit.s2;
287    trace.final_noise1 = unit.noise1;
288    trace.final_noise2 = unit.noise2;
289    Ok(trace)
290}
291
292#[cfg(test)]
293mod tests {
294    use super::*;
295
296    #[test]
297    fn one_step_matches_appendix_euler_and_ou_equations() {
298        let mut unit = WongWangUnit::with_parameters(
299            0.1, 0.2, 0.01, -0.02, 0.1, 0.002, 0.641, 0.2609, 0.0497, 0.3255, 0.02, 0.0001, 3,
300        )
301        .unwrap();
302        let old = (unit.s1, unit.s2, unit.noise1, unit.noise2);
303        let (rate1, rate2) = unit
304            .step_with_gaussian_samples(0.03, -0.01, 0.5, -1.0)
305            .unwrap();
306        let expected_s1 = old.0 + 0.0001 * (-old.0 / 0.1 + (1.0 - old.0) * 0.641 * rate1);
307        let expected_s2 = old.1 + 0.0001 * (-old.1 / 0.1 + (1.0 - old.1) * 0.641 * rate2);
308        let scale = (0.0001_f64 / 0.002).sqrt() * 0.02;
309        assert_eq!(unit.s1, expected_s1);
310        assert_eq!(unit.s2, expected_s2);
311        assert_eq!(unit.noise1, old.2 - 0.05 * old.2 + scale * 0.5);
312        assert_eq!(unit.noise2, old.3 - 0.05 * old.3 - scale);
313    }
314
315    #[test]
316    fn batch_matches_scalar_and_preserves_empty_initial_state() {
317        let empty = simulate(
318            0.2,
319            0.3,
320            0.01,
321            -0.02,
322            0.1,
323            0.002,
324            0.641,
325            0.2609,
326            0.0497,
327            0.3255,
328            0.02,
329            0.0001,
330            &[],
331            &[],
332            &[],
333        )
334        .unwrap();
335        assert!(empty.s1.is_empty());
336        assert_eq!(empty.final_s1, 0.2);
337        assert_eq!(empty.final_noise2, -0.02);
338
339        let stim1 = [0.02, 0.01, -0.01];
340        let stim2 = [-0.01, 0.0, 0.03];
341        let xi = [0.5, -1.0, 0.25, 0.75, -0.5, 0.0];
342        let batch = simulate(
343            0.2, 0.3, 0.01, -0.02, 0.12, 0.003, 0.7, 0.28, 0.06, 0.31, 0.015, 0.0002, &stim1,
344            &stim2, &xi,
345        )
346        .unwrap();
347        let mut scalar = WongWangUnit::with_parameters(
348            0.2, 0.3, 0.01, -0.02, 0.12, 0.003, 0.7, 0.28, 0.06, 0.31, 0.015, 0.0002, 0,
349        )
350        .unwrap();
351        for step in 0..stim1.len() {
352            let rates = scalar
353                .step_with_gaussian_samples(
354                    stim1[step],
355                    stim2[step],
356                    xi[2 * step],
357                    xi[2 * step + 1],
358                )
359                .unwrap();
360            assert_eq!(batch.s1[step], scalar.s1);
361            assert_eq!(batch.s2[step], scalar.s2);
362            assert_eq!(batch.noise1[step], scalar.noise1);
363            assert_eq!(batch.noise2[step], scalar.noise2);
364            assert_eq!((batch.rate1[step], batch.rate2[step]), rates);
365        }
366        assert_eq!(batch.final_s1, scalar.s1);
367        assert_eq!(batch.final_s2, scalar.s2);
368        assert_eq!(batch.final_noise1, scalar.noise1);
369        assert_eq!(batch.final_noise2, scalar.noise2);
370    }
371
372    #[test]
373    fn invalid_input_is_rejected_before_state_commit() {
374        let mut unit = WongWangUnit::new(5);
375        let before = (unit.s1, unit.s2, unit.noise1, unit.noise2);
376        assert!(unit
377            .step_with_gaussian_samples(f64::NAN, 0.0, 0.0, 0.0)
378            .is_err());
379        assert_eq!((unit.s1, unit.s2, unit.noise1, unit.noise2), before);
380        assert!(simulate(
381            0.1,
382            0.1,
383            0.0,
384            0.0,
385            0.1,
386            0.002,
387            0.641,
388            0.2609,
389            0.0497,
390            0.3255,
391            0.02,
392            0.0001,
393            &[0.0],
394            &[0.0],
395            &[],
396        )
397        .is_err());
398    }
399
400    #[test]
401    fn reset_preserves_configured_parameters() {
402        let mut unit = WongWangUnit::with_parameters(
403            0.2, 0.3, 0.01, -0.02, 0.12, 0.003, 0.7, 0.3, 0.04, 0.31, 0.03, 0.0002, 7,
404        )
405        .unwrap();
406        unit.reset();
407        assert_eq!(
408            (unit.s1, unit.s2, unit.noise1, unit.noise2),
409            (0.1, 0.1, 0.0, 0.0)
410        );
411        assert_eq!((unit.tau_s, unit.tau_ampa, unit.dt), (0.12, 0.003, 0.0002));
412    }
413}
sc_neurocore_engine/wong_wang.rs

sc_neurocore_engine/
wong_wang.rs
