sc_neurocore_engine/
ei_network.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 — Floating-point E-I balanced network for Studio
8
9//! Fused E-I balanced LIF network simulation in f64 arithmetic.
10//!
11//! Runs a complete excitatory-inhibitory network with CSR weight matrix,
12//! Poisson external drive, and spike-scatter coupling in a single Rust
13//! call — no per-step Python overhead.
14//!
15//! Synapses are delta synapses in the sense of Brunel (2000, J. Comput.
16//! Neurosci. 8:183): a presynaptic spike moves the postsynaptic membrane by
17//! its weight in millivolts, one step later. Every neuron also receives
18//! [`EXTERNAL_SYNAPSES`] independent excitatory Poisson inputs, each firing at
19//! `ext_rate` Hz and moving the membrane by [`EXTERNAL_WEIGHT_MV`]. The
20//! external drive alone reaches threshold at about 9.4 Hz per input
21//! (15 mV / (0.1 mV x 800 x 20 ms)).
22//!
23//! An earlier version multiplied each jump by the step size, so a spike moved
24//! the membrane by `weight x dt` (0.01 mV for the default E-to-E weight at
25//! dt = 0.1 ms), and the external drive was one input per neuron: no setting
26//! of the Studio's controls produced a single spike.
27
28use rand::{RngExt, SeedableRng};
29use rand_distr::{Distribution, Poisson};
30use rand_xoshiro::Xoshiro256PlusPlus;
31
32/// Independent external excitatory Poisson inputs per neuron.
33pub const EXTERNAL_SYNAPSES: f64 = 800.0;
34
35/// Membrane jump caused by one external input spike, in millivolts.
36pub const EXTERNAL_WEIGHT_MV: f64 = 0.1;
37
38/// Result of an E-I network simulation.
39pub struct EIResult {
40    pub spike_times: Vec<f64>,
41    pub spike_neurons: Vec<u32>,
42    pub n_exc: u32,
43    pub n_inh: u32,
44    pub rate_time: Vec<f64>,
45    pub exc_rates: Vec<f64>,
46    pub inh_rates: Vec<f64>,
47    pub mean_exc_rate: f64,
48    pub mean_inh_rate: f64,
49}
50
51/// Run a complete E-I balanced LIF network simulation.
52///
53/// All computation happens in Rust — connectivity build, Poisson input,
54/// Euler integration, spike detection, rate binning.
55#[allow(clippy::too_many_arguments)]
56pub fn simulate_ei(
57    n_exc: usize,
58    n_inh: usize,
59    w_ee: f64,
60    w_ei: f64,
61    w_ie: f64,
62    w_ii: f64,
63    p_conn: f64,
64    ext_rate: f64,
65    duration: f64,
66    dt: f64,
67    seed: u64,
68) -> EIResult {
69    let n = n_exc + n_inh;
70    let n_steps = ((duration / dt) as usize).min(50_000);
71    let mut rng = Xoshiro256PlusPlus::seed_from_u64(seed);
72
73    // LIF params
74    let tau_m = 20.0_f64;
75    let v_rest = -65.0_f64;
76    let v_threshold = -50.0_f64;
77    let v_reset = -65.0_f64;
78    let tau_ref = 2.0_f64;
79
80    // Build CSR weight matrix
81    let mut row_offsets = Vec::with_capacity(n + 1);
82    let mut col_indices = Vec::new();
83    let mut values = Vec::new();
84    row_offsets.push(0usize);
85
86    for i in 0..n {
87        let i_exc = i < n_exc;
88        for j in 0..n {
89            if i == j {
90                continue;
91            }
92            if rng.random::<f64>() >= p_conn {
93                continue;
94            }
95            let j_exc = j < n_exc;
96            let w = match (i_exc, j_exc) {
97                (true, true) => w_ee,
98                (true, false) => w_ie,
99                (false, true) => -w_ei,
100                (false, false) => -w_ii,
101            };
102            col_indices.push(j);
103            values.push(w);
104        }
105        row_offsets.push(col_indices.len());
106    }
107
108    // State arrays
109    let mut v = vec![v_rest; n];
110    let mut refractory = vec![0.0_f64; n];
111    let mut prev_spiked = vec![false; n];
112
113    // Recording
114    let mut spike_times: Vec<f64> = Vec::new();
115    let mut spike_neurons: Vec<u32> = Vec::new();
116    let bin_size = (n_steps / 100).max(1);
117    let n_bins = n_steps / bin_size;
118    let mut exc_rates = vec![0.0_f64; n_bins];
119    let mut inh_rates = vec![0.0_f64; n_bins];
120    let mut exc_bin = 0u32;
121    let mut inh_bin = 0u32;
122
123    // Expected external input spikes per neuron per step. `Poisson` samples
124    // exactly for any mean; the multiply-uniforms loop it replaces never
125    // terminated once exp(-mean) underflowed to zero.
126    let ext_lambda = EXTERNAL_SYNAPSES * ext_rate.max(0.0) * dt / 1000.0;
127    let ext_inputs = if ext_lambda > 0.0 {
128        Poisson::new(ext_lambda).ok()
129    } else {
130        None
131    };
132    let mut exc_spikes = 0usize;
133    let mut inh_spikes = 0usize;
134
135    for t in 0..n_steps {
136        // Decay refractory
137        for r in refractory.iter_mut() {
138            *r = (*r - dt).max(0.0);
139        }
140
141        // Synaptic input from previous spikes (CSR scatter)
142        let mut syn = vec![0.0_f64; n];
143        for i in 0..n {
144            if !prev_spiked[i] {
145                continue;
146            }
147            let start = row_offsets[i];
148            let end = row_offsets[i + 1];
149            for k in start..end {
150                syn[col_indices[k]] += values[k];
151            }
152        }
153
154        // External Poisson + Euler step
155        for i in 0..n {
156            if refractory[i] > 0.0 {
157                continue;
158            }
159            // Knuth Poisson sampling (fast for small lambda)
160            let external = match &ext_inputs {
161                Some(poisson) => poisson.sample(&mut rng) * EXTERNAL_WEIGHT_MV,
162                None => 0.0,
163            };
164            v[i] += -(v[i] - v_rest) / tau_m * dt + external + syn[i];
165        }
166
167        // Spike detection
168        prev_spiked.fill(false);
169        for i in 0..n {
170            if refractory[i] > 0.0 {
171                continue;
172            }
173            if v[i] >= v_threshold {
174                v[i] = v_reset;
175                refractory[i] = tau_ref;
176                prev_spiked[i] = true;
177                spike_times.push(t as f64 * dt);
178                spike_neurons.push(i as u32);
179                if i < n_exc {
180                    exc_bin += 1;
181                    exc_spikes += 1;
182                } else {
183                    inh_bin += 1;
184                    inh_spikes += 1;
185                }
186            }
187        }
188
189        // Rate binning
190        if (t + 1) % bin_size == 0 {
191            let bi = t / bin_size;
192            if bi < n_bins {
193                let bin_t = bin_size as f64 * dt / 1000.0;
194                exc_rates[bi] = exc_bin as f64 / n_exc.max(1) as f64 / bin_t.max(0.001);
195                inh_rates[bi] = inh_bin as f64 / n_inh.max(1) as f64 / bin_t.max(0.001);
196            }
197            exc_bin = 0;
198            inh_bin = 0;
199        }
200    }
201
202    let rate_time: Vec<f64> = (0..n_bins)
203        .map(|i| i as f64 * bin_size as f64 * dt)
204        .collect();
205
206    // Mean rate per neuron over the simulated time. It used to average only
207    // the rate bins that held a spike, which overstated sparse activity.
208    let simulated_s = (n_steps as f64 * dt / 1000.0).max(f64::MIN_POSITIVE);
209    let mean_exc = exc_spikes as f64 / n_exc.max(1) as f64 / simulated_s;
210    let mean_inh = inh_spikes as f64 / n_inh.max(1) as f64 / simulated_s;
211
212    EIResult {
213        spike_times,
214        spike_neurons,
215        n_exc: n_exc as u32,
216        n_inh: n_inh as u32,
217        rate_time,
218        exc_rates,
219        inh_rates,
220        mean_exc_rate: (mean_exc * 10.0).round() / 10.0,
221        mean_inh_rate: (mean_inh * 10.0).round() / 10.0,
222    }
223}
224
225#[cfg(test)]
226mod tests {
227    use super::*;
228
229    const DEFAULT: (f64, f64, f64, f64, f64) = (0.1, 0.4, 0.1, 0.4, 0.2);
230
231    fn run(ext_rate: f64, duration: f64, dt: f64) -> EIResult {
232        let (w_ee, w_ei, w_ie, w_ii, p_conn) = DEFAULT;
233        simulate_ei(
234            80, 20, w_ee, w_ei, w_ie, w_ii, p_conn, ext_rate, duration, dt, 42,
235        )
236    }
237
238    #[test]
239    fn ei_network_runs_without_panic() {
240        let r = simulate_ei(20, 5, 0.1, 0.4, 0.1, 0.4, 0.2, 10.0, 50.0, 0.1, 42);
241        assert_eq!(r.n_exc, 20);
242        assert_eq!(r.n_inh, 5);
243        assert!(!r.rate_time.is_empty());
244        assert_eq!(r.exc_rates.len(), r.rate_time.len());
245    }
246
247    #[test]
248    fn default_drive_above_threshold_fires() {
249        // 12 Hz per input is 1.3x the 9.4 Hz threshold rate of the drive.
250        let r = run(12.0, 500.0, 0.1);
251        assert!(
252            r.mean_exc_rate > 5.0 && r.mean_exc_rate < 100.0,
253            "{}",
254            r.mean_exc_rate
255        );
256        assert!(
257            r.mean_inh_rate > 5.0 && r.mean_inh_rate < 100.0,
258            "{}",
259            r.mean_inh_rate
260        );
261    }
262
263    #[test]
264    fn drive_well_below_threshold_stays_silent() {
265        // Mean drive 3.2 mV against a 15 mV gap, fluctuations of about 0.3 mV.
266        assert!(run(2.0, 500.0, 0.1).spike_times.is_empty());
267    }
268
269    #[test]
270    fn rate_does_not_depend_on_the_step_size() {
271        // A spike moves the membrane by its weight, not by weight x dt.
272        let coarse = run(12.0, 1000.0, 0.1).mean_exc_rate;
273        let fine = run(12.0, 1000.0, 0.05).mean_exc_rate;
274        assert!(
275            (coarse - fine).abs() < 0.2 * coarse.max(fine),
276            "{coarse} vs {fine}"
277        );
278    }
279
280    #[test]
281    fn mean_rate_counts_silent_time() {
282        let r = run(12.0, 500.0, 0.1);
283        let spikes = r.spike_neurons.iter().filter(|&&i| i < 80).count() as f64;
284        assert!((r.mean_exc_rate - spikes / 80.0 / 0.5).abs() < 0.051);
285    }
286
287    #[test]
288    fn very_large_external_mean_terminates() {
289        // 800 inputs x 100 Hz x 5 ms = 400 expected events per step.
290        let r = run(100.0, 50.0, 5.0);
291        assert!(!r.spike_times.is_empty());
292    }
293}
sc_neurocore_engine/ei_network.rs

sc_neurocore_engine/
ei_network.rs
