sc_neurocore_engine/neurons/
multi_compartment.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 — Multi-compartment neuron models
8
9//! Multi-compartment neuron models.
10
11#[allow(unused_imports)]
12use super::biophysical::safe_rate;
13
14/// Pinsky-Rinzel 1994 — 2-compartment pyramidal cell.
15#[derive(Clone, Debug)]
16pub struct PinskyRinzelNeuron {
17    pub v_s: f64,
18    pub v_d: f64,
19    pub h: f64,
20    pub n: f64,
21    pub s: f64,
22    pub c: f64,
23    pub q: f64,
24    pub gc: f64,
25    pub p: f64,
26    pub g_na: f64,
27    pub g_kdr: f64,
28    pub g_ca: f64,
29    pub g_kahp: f64,
30    pub g_l: f64,
31    pub e_na: f64,
32    pub e_k: f64,
33    pub e_ca: f64,
34    pub e_l: f64,
35    pub dt: f64,
36    pub v_threshold: f64,
37}
38
39impl PinskyRinzelNeuron {
40    pub fn new() -> Self {
41        Self {
42            v_s: -60.0,
43            v_d: -60.0,
44            h: 0.9,
45            n: 0.1,
46            s: 0.0,
47            c: 0.0,
48            q: 0.0,
49            gc: 2.1,
50            p: 0.5,
51            g_na: 30.0,
52            g_kdr: 15.0,
53            g_ca: 10.0,
54            g_kahp: 0.8,
55            g_l: 0.1,
56            e_na: 60.0,
57            e_k: -75.0,
58            e_ca: 80.0,
59            e_l: -60.0,
60            dt: 0.02,
61            v_threshold: -20.0,
62        }
63    }
64    pub fn step(&mut self, current_soma: f64, current_dend: f64) -> i32 {
65        let v_prev = self.v_s;
66        let am = safe_rate(0.32, 54.0, self.v_s, 4.0, 8.0);
67        let bm = safe_rate(-0.28, 27.0, self.v_s, -5.0, 5.6);
68        let m_inf = am / (am + bm);
69        let ah = 0.128 * (-(self.v_s + 50.0) / 18.0).exp();
70        let bh = 4.0 / (1.0 + (-(self.v_s + 27.0) / 5.0).exp());
71        let an = safe_rate(0.032, 52.0, self.v_s, 5.0, 0.32);
72        let bn = 0.5 * (-(self.v_s + 57.0) / 40.0).exp();
73        let s_inf = 1.0 / (1.0 + (-(self.v_d + 20.0) / 9.0).exp());
74        let i_na = self.g_na * m_inf.powi(2) * self.h * (self.v_s - self.e_na);
75        let i_kdr = self.g_kdr * self.n.powi(2) * (self.v_s - self.e_k);
76        let i_ls = self.g_l * (self.v_s - self.e_l);
77        let i_ds = (self.gc / self.p) * (self.v_s - self.v_d);
78        let i_ca = self.g_ca * self.s.powi(2) * (self.v_d - self.e_ca);
79        let i_kahp = self.g_kahp * self.q * (self.v_d - self.e_k);
80        let i_ld = self.g_l * (self.v_d - self.e_l);
81        let i_sd = (self.gc / (1.0 - self.p)) * (self.v_d - self.v_s);
82        self.v_s += (-i_na - i_kdr - i_ls - i_ds + current_soma / self.p) * self.dt;
83        self.v_d += (-i_ca - i_kahp - i_ld - i_sd + current_dend / (1.0 - self.p)) * self.dt;
84        self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
85        self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
86        self.s += ((s_inf - self.s) / 5.0) * self.dt;
87        self.c = (self.c + (-0.13 * i_ca - 0.075 * self.c) * self.dt).max(0.0);
88        let q_inf = (self.c / (self.c + 2.0)).min(1.0);
89        self.q += ((q_inf - self.q) / 100.0) * self.dt;
90        if self.v_s >= self.v_threshold && v_prev < self.v_threshold {
91            1
92        } else {
93            0
94        }
95    }
96    pub fn reset(&mut self) {
97        self.v_s = -60.0;
98        self.v_d = -60.0;
99        self.h = 0.9;
100        self.n = 0.1;
101        self.s = 0.0;
102        self.c = 0.0;
103        self.q = 0.0;
104    }
105}
106impl Default for PinskyRinzelNeuron {
107    fn default() -> Self {
108        Self::new()
109    }
110}
111
112/// Hay et al. 2011 — Layer 5 thick-tufted pyramidal (3-compartment reduced).
113#[derive(Clone, Debug)]
114pub struct HayL5PyramidalNeuron {
115    pub v_s: f64,
116    pub v_t: f64,
117    pub v_a: f64,
118    pub h_na: f64,
119    pub n_k: f64,
120    pub m_ca: f64,
121    pub h_ca: f64,
122    pub m_ih: f64,
123    pub ca_a: f64,
124    pub g_na: f64,
125    pub g_k: f64,
126    pub g_l_s: f64,
127    pub g_ca_t: f64,
128    pub g_ih: f64,
129    pub g_l_t: f64,
130    pub g_ca_a: f64,
131    pub g_kca: f64,
132    pub g_l_a: f64,
133    pub g_st: f64,
134    pub g_ta: f64,
135    pub p_s: f64,
136    pub p_t: f64,
137    pub p_a: f64,
138    pub e_na: f64,
139    pub e_k: f64,
140    pub e_ca: f64,
141    pub e_ih: f64,
142    pub e_l: f64,
143    pub ca_decay: f64,
144    pub f_ca: f64,
145    pub c_m: f64,
146    pub dt: f64,
147    pub v_threshold: f64,
148}
149
150impl HayL5PyramidalNeuron {
151    pub fn new() -> Self {
152        Self {
153            v_s: -75.0,
154            v_t: -75.0,
155            v_a: -75.0,
156            h_na: 0.9,
157            n_k: 0.1,
158            m_ca: 0.0,
159            h_ca: 1.0,
160            m_ih: 0.0,
161            ca_a: 0.0001,
162            g_na: 300.0,
163            g_k: 40.0,
164            g_l_s: 0.03,
165            g_ca_t: 2.0,
166            g_ih: 0.02,
167            g_l_t: 0.03,
168            g_ca_a: 1.5,
169            g_kca: 2.5,
170            g_l_a: 0.03,
171            g_st: 1.5,
172            g_ta: 0.8,
173            p_s: 0.15,
174            p_t: 0.25,
175            p_a: 0.60,
176            e_na: 50.0,
177            e_k: -85.0,
178            e_ca: 140.0,
179            e_ih: -45.0,
180            e_l: -75.0,
181            ca_decay: 200.0,
182            f_ca: 0.0002,
183            c_m: 1.0,
184            dt: 0.025,
185            v_threshold: -30.0,
186        }
187    }
188    pub fn step(&mut self, current_soma: f64, current_tuft: f64) -> i32 {
189        let v_s_prev = self.v_s;
190        for _ in 0..4 {
191            let m_na_inf = 1.0 / (1.0 + (-(self.v_s + 38.0) / 7.0).exp());
192            let h_na_inf = 1.0 / (1.0 + ((self.v_s + 65.0) / 6.0).exp());
193            let n_k_inf = 1.0 / (1.0 + (-(self.v_s + 25.0) / 12.0).exp());
194            let tau_h = 0.5 + 14.0 / (1.0 + ((self.v_s + 35.0) / 10.0).exp());
195            let tau_n = 1.0 + 5.0 / (1.0 + ((self.v_s + 30.0) / 10.0).exp());
196            self.h_na += (h_na_inf - self.h_na) / tau_h * self.dt;
197            self.n_k += (n_k_inf - self.n_k) / tau_n * self.dt;
198            let i_na = self.g_na * m_na_inf.powi(3) * self.h_na * (self.v_s - self.e_na);
199            let i_k = self.g_k * self.n_k.powi(4) * (self.v_s - self.e_k);
200            let i_l_s = self.g_l_s * (self.v_s - self.e_l);
201            let i_st = self.g_st * (self.v_s - self.v_t) / self.p_s;
202            let m_ca_inf = 1.0 / (1.0 + (-(self.v_t + 27.0) / 7.0).exp());
203            let h_ca_inf = 1.0 / (1.0 + ((self.v_t + 52.0) / 5.0).exp());
204            let m_ih_inf = 1.0 / (1.0 + ((self.v_t + 75.0) / 5.5).exp());
205            self.m_ca += (m_ca_inf - self.m_ca) / 1.0 * self.dt;
206            self.h_ca += (h_ca_inf - self.h_ca) / 20.0 * self.dt;
207            self.m_ih += (m_ih_inf - self.m_ih) / 50.0 * self.dt;
208            let i_ca_t = self.g_ca_t * self.m_ca.powi(2) * self.h_ca * (self.v_t - self.e_ca);
209            let i_ih = self.g_ih * self.m_ih * (self.v_t - self.e_ih);
210            let i_l_t = self.g_l_t * (self.v_t - self.e_l);
211            let i_ts = self.g_st * (self.v_t - self.v_s) / self.p_t;
212            let i_ta = self.g_ta * (self.v_t - self.v_a) / self.p_t;
213            let m_ca_a_inf = 1.0 / (1.0 + (-(self.v_a + 30.0) / 5.0).exp());
214            let kca_act = self.ca_a / (self.ca_a + 0.001);
215            let i_ca_a = self.g_ca_a * m_ca_a_inf.powi(2) * (self.v_a - self.e_ca);
216            let i_kca = self.g_kca * kca_act * (self.v_a - self.e_k);
217            let i_l_a = self.g_l_a * (self.v_a - self.e_l);
218            let i_at = self.g_ta * (self.v_a - self.v_t) / self.p_a;
219            self.v_s += (-i_na - i_k - i_l_s - i_st + current_soma / self.p_s) / self.c_m * self.dt;
220            self.v_t += (-i_ca_t - i_ih - i_l_t - i_ts - i_ta) / self.c_m * self.dt;
221            self.v_a +=
222                (-i_ca_a - i_kca - i_l_a - i_at + current_tuft / self.p_a) / self.c_m * self.dt;
223            self.ca_a =
224                (self.ca_a + (-self.f_ca * i_ca_a - self.ca_a / self.ca_decay) * self.dt).max(0.0);
225        }
226        if self.v_s >= self.v_threshold && v_s_prev < self.v_threshold {
227            1
228        } else {
229            0
230        }
231    }
232    pub fn reset(&mut self) {
233        self.v_s = -75.0;
234        self.v_t = -75.0;
235        self.v_a = -75.0;
236        self.h_na = 0.9;
237        self.n_k = 0.1;
238        self.m_ca = 0.0;
239        self.h_ca = 1.0;
240        self.m_ih = 0.0;
241        self.ca_a = 0.0001;
242    }
243}
244impl Default for HayL5PyramidalNeuron {
245    fn default() -> Self {
246        Self::new()
247    }
248}
249
250/// Marder STG — stomatogastric ganglion LP neuron with 7 currents.
251#[derive(Clone, Debug)]
252pub struct MarderSTGNeuron {
253    pub v: f64,
254    pub m_na: f64,
255    pub h_na: f64,
256    pub m_cat: f64,
257    pub h_cat: f64,
258    pub m_cas: f64,
259    pub m_a: f64,
260    pub h_a: f64,
261    pub m_kd: f64,
262    pub m_h: f64,
263    pub ca: f64,
264    pub g_na: f64,
265    pub g_cat: f64,
266    pub g_cas: f64,
267    pub g_a: f64,
268    pub g_kd: f64,
269    pub g_kca: f64,
270    pub g_h: f64,
271    pub g_l: f64,
272    pub e_na: f64,
273    pub e_k: f64,
274    pub e_ca: f64,
275    pub e_h: f64,
276    pub e_l: f64,
277    pub dt: f64,
278    pub v_threshold: f64,
279}
280
281impl MarderSTGNeuron {
282    pub fn new() -> Self {
283        Self {
284            v: -60.0,
285            m_na: 0.0,
286            h_na: 0.9,
287            m_cat: 0.0,
288            h_cat: 0.9,
289            m_cas: 0.0,
290            m_a: 0.0,
291            h_a: 0.9,
292            m_kd: 0.0,
293            m_h: 0.0,
294            ca: 0.05,
295            g_na: 400.0,
296            g_cat: 2.5,
297            g_cas: 6.0,
298            g_a: 50.0,
299            g_kd: 100.0,
300            g_kca: 25.0,
301            g_h: 0.02,
302            g_l: 0.01,
303            e_na: 50.0,
304            e_k: -80.0,
305            e_ca: 120.0,
306            e_h: -20.0,
307            e_l: -50.0,
308            dt: 0.05,
309            v_threshold: -20.0,
310        }
311    }
312    pub fn step(&mut self, current: f64) -> i32 {
313        let v_prev = self.v;
314        let b = |v: f64, vh: f64, s: f64| 1.0 / (1.0 + (-(v - vh) / s).exp());
315        self.m_na += (b(self.v, -25.5, 5.29) - self.m_na) / 1.32 * self.dt;
316        self.h_na += (b(self.v, -48.9, -5.18) - self.h_na)
317            / (0.67 * (1.0 + ((self.v + 62.9) / -10.0).exp()) + 1.5)
318            * self.dt;
319        self.m_cat += (b(self.v, -27.1, 7.2) - self.m_cat) / 21.7 * self.dt;
320        self.h_cat += (b(self.v, -32.1, -5.5) - self.h_cat) / 105.0 * self.dt;
321        self.m_cas += (b(self.v, -33.0, 8.1) - self.m_cas) / 14.0 * self.dt;
322        self.m_a += (b(self.v, -27.2, 8.7) - self.m_a) / 11.6 * self.dt;
323        self.h_a += (b(self.v, -56.9, -4.9) - self.h_a) / 38.6 * self.dt;
324        self.m_kd += (b(self.v, -12.3, 11.8) - self.m_kd) / 7.2 * self.dt;
325        self.m_h += (b(self.v, -70.0, -6.0) - self.m_h) / 272.0 * self.dt;
326        let kca_act = (self.ca / (self.ca + 3.0)).min(1.0);
327        let i_na = self.g_na * self.m_na.powi(3) * self.h_na * (self.v - self.e_na);
328        let i_cat = self.g_cat * self.m_cat.powi(3) * self.h_cat * (self.v - self.e_ca);
329        let i_cas = self.g_cas * self.m_cas.powi(3) * (self.v - self.e_ca);
330        let i_a = self.g_a * self.m_a.powi(3) * self.h_a * (self.v - self.e_k);
331        let i_kd = self.g_kd * self.m_kd.powi(4) * (self.v - self.e_k);
332        let i_kca = self.g_kca * kca_act.powi(4) * (self.v - self.e_k);
333        let i_h = self.g_h * self.m_h * (self.v - self.e_h);
334        let i_l = self.g_l * (self.v - self.e_l);
335        self.v += (-i_na - i_cat - i_cas - i_a - i_kd - i_kca - i_h - i_l + current) * self.dt;
336        let i_ca_total = i_cat + i_cas;
337        self.ca = (self.ca + (-0.0001 * i_ca_total - 0.01 * self.ca) * self.dt).max(0.0);
338        if self.v >= self.v_threshold && v_prev < self.v_threshold {
339            1
340        } else {
341            0
342        }
343    }
344    pub fn reset(&mut self) {
345        self.v = -60.0;
346        self.m_na = 0.0;
347        self.h_na = 0.9;
348        self.m_cat = 0.0;
349        self.h_cat = 0.9;
350        self.m_cas = 0.0;
351        self.m_a = 0.0;
352        self.h_a = 0.9;
353        self.m_kd = 0.0;
354        self.m_h = 0.0;
355        self.ca = 0.05;
356    }
357}
358impl Default for MarderSTGNeuron {
359    fn default() -> Self {
360        Self::new()
361    }
362}
363
364/// Rall cable — N-compartment passive dendrite model. Rall 1964.
365#[derive(Clone, Debug)]
366pub struct RallCableNeuron {
367    pub v: Vec<f64>,
368    pub n_comp: usize,
369    pub tau_m: f64,
370    pub v_rest: f64,
371    pub g_ratio: f64,
372    pub v_threshold: f64,
373    pub v_reset: f64,
374    pub dt: f64,
375}
376
377impl RallCableNeuron {
378    pub fn new(n_comp: usize) -> Self {
379        Self {
380            v: vec![-65.0; n_comp],
381            n_comp,
382            tau_m: 20.0,
383            v_rest: -65.0,
384            g_ratio: 0.5,
385            v_threshold: -50.0,
386            v_reset: -65.0,
387            dt: 0.1,
388        }
389    }
390    pub fn step(&mut self, current: f64) -> i32 {
391        let n = self.n_comp;
392        let mut dv = vec![0.0; n];
393        for i in 0..n {
394            let leak = -(self.v[i] - self.v_rest) / self.tau_m;
395            let left = if i > 0 {
396                self.g_ratio * (self.v[i - 1] - self.v[i])
397            } else {
398                0.0
399            };
400            let right = if i + 1 < n {
401                self.g_ratio * (self.v[i + 1] - self.v[i])
402            } else {
403                0.0
404            };
405            let inj = if i == n - 1 { current } else { 0.0 };
406            dv[i] = (leak + left + right + inj) * self.dt;
407        }
408        for i in 0..n {
409            self.v[i] += dv[i];
410        }
411        if self.v[0] >= self.v_threshold {
412            self.v[0] = self.v_reset;
413            1
414        } else {
415            0
416        }
417    }
418    pub fn reset(&mut self) {
419        self.v.fill(self.v_rest);
420    }
421}
422
423/// Booth-Rinzel — 2-compartment motoneuron with bistability. Booth et al. 1997.
424#[derive(Clone, Debug)]
425pub struct BoothRinzelNeuron {
426    pub vs: f64,
427    pub vd: f64,
428    pub h: f64,
429    pub n: f64,
430    pub q: f64,
431    pub ca: f64,
432    pub p: f64,
433    pub gc: f64,
434    pub g_na: f64,
435    pub g_k: f64,
436    pub g_ca: f64,
437    pub g_kca: f64,
438    pub g_l: f64,
439    pub e_na: f64,
440    pub e_k: f64,
441    pub e_ca: f64,
442    pub e_l: f64,
443    pub dt: f64,
444    pub v_threshold: f64,
445}
446
447impl BoothRinzelNeuron {
448    pub fn new() -> Self {
449        Self {
450            vs: -65.0,
451            vd: -65.0,
452            h: 0.9,
453            n: 0.0,
454            q: 0.0,
455            ca: 0.0,
456            p: 0.5,
457            gc: 0.1,
458            g_na: 120.0,
459            g_k: 20.0,
460            g_ca: 14.0,
461            g_kca: 5.0,
462            g_l: 0.51,
463            e_na: 55.0,
464            e_k: -80.0,
465            e_ca: 80.0,
466            e_l: -60.0,
467            dt: 0.025,
468            v_threshold: -20.0,
469        }
470    }
471    pub fn step(&mut self, current: f64) -> i32 {
472        let v_prev = self.vs;
473        for _ in 0..4 {
474            let m_inf = 1.0 / (1.0 + (-(self.vs + 35.0) / 7.8).exp());
475            let h_inf = 1.0 / (1.0 + ((self.vs + 55.0) / 7.0).exp());
476            let n_inf = 1.0 / (1.0 + (-(self.vs + 28.0) / 15.0).exp());
477            let s_inf = 1.0 / (1.0 + (-(self.vd + 22.0) / 5.0).exp());
478            let q_inf = 1.0 / (1.0 + (-(self.vd + 35.0) / 2.0).exp());
479            let tau_h = (30.0
480                / (((self.vs + 50.0) / 15.0).exp() + ((-(self.vs + 50.0)) / 16.0).exp() + 1e-12))
481                .max(0.01);
482            let tau_n = (7.0
483                / (((self.vs + 40.0) / 40.0).exp() + ((-(self.vs + 40.0)) / 50.0).exp() + 1e-12))
484                .max(0.01);
485            self.h = (self.h + (h_inf - self.h) / tau_h * self.dt).clamp(0.0, 1.0);
486            self.n = (self.n + (n_inf - self.n) / tau_n * self.dt).clamp(0.0, 1.0);
487            self.q = (self.q + (q_inf - self.q) / 400.0 * self.dt).clamp(0.0, 1.0);
488            let chi = (self.ca / 250.0).min(1.0);
489            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.vs - self.e_na);
490            let i_k = self.g_k * self.n.powi(4) * (self.vs - self.e_k);
491            let i_ls = self.g_l * (self.vs - self.e_l);
492            let i_sd = (self.gc / self.p) * (self.vs - self.vd);
493            let i_ca = self.g_ca * s_inf.powi(2) * (self.vd - self.e_ca);
494            let i_kca = self.g_kca * chi * (self.vd - self.e_k);
495            let i_ld = self.g_l * (self.vd - self.e_l);
496            let i_ds = (self.gc / (1.0 - self.p)) * (self.vd - self.vs);
497            self.vs = (self.vs + (-i_na - i_k - i_ls - i_sd + current / self.p) * self.dt)
498                .clamp(-200.0, 100.0);
499            self.vd = (self.vd + (-i_ca - i_kca - i_ld - i_ds) * self.dt).clamp(-200.0, 100.0);
500            self.ca = (self.ca + (0.0025 * (-0.009 * i_ca) - 0.18 * self.ca) * self.dt).max(0.0);
501        }
502        if self.vs >= self.v_threshold && v_prev < self.v_threshold {
503            1
504        } else {
505            0
506        }
507    }
508    pub fn reset(&mut self) {
509        self.vs = -65.0;
510        self.vd = -65.0;
511        self.h = 0.9;
512        self.n = 0.0;
513        self.q = 0.0;
514        self.ca = 0.0;
515    }
516}
517impl Default for BoothRinzelNeuron {
518    fn default() -> Self {
519        Self::new()
520    }
521}
522
523/// Dendrify — two-compartment with active dendritic spike (NMDA-like plateau).
524#[derive(Clone, Debug)]
525pub struct DendrifyNeuron {
526    pub v_s: f64,
527    pub v_d: f64,
528    pub d_active: bool,
529    pub d_timer: f64,
530    pub tau_s: f64,
531    pub tau_d: f64,
532    pub g_c: f64,
533    pub d_threshold: f64,
534    pub d_amplitude: f64,
535    pub d_duration: f64,
536    pub v_rest: f64,
537    pub v_threshold: f64,
538    pub v_reset: f64,
539    pub dt: f64,
540}
541
542impl DendrifyNeuron {
543    pub fn new() -> Self {
544        Self {
545            v_s: -65.0,
546            v_d: -65.0,
547            d_active: false,
548            d_timer: 0.0,
549            tau_s: 10.0,
550            tau_d: 20.0,
551            g_c: 0.8,
552            d_threshold: -35.0,
553            d_amplitude: 30.0,
554            d_duration: 10.0,
555            v_rest: -65.0,
556            v_threshold: -50.0,
557            v_reset: -65.0,
558            dt: 0.1,
559        }
560    }
561    pub fn step(&mut self, current: f64) -> i32 {
562        let d_input = if self.d_active { self.d_amplitude } else { 0.0 };
563        self.v_d += (-(self.v_d - self.v_rest) + current + d_input
564            - self.g_c * (self.v_d - self.v_s))
565            / self.tau_d
566            * self.dt;
567        self.v_s +=
568            (-(self.v_s - self.v_rest) + self.g_c * (self.v_d - self.v_s)) / self.tau_s * self.dt;
569        if self.d_active {
570            self.d_timer -= self.dt;
571            if self.d_timer <= 0.0 {
572                self.d_active = false;
573            }
574        } else if self.v_d >= self.d_threshold {
575            self.d_active = true;
576            self.d_timer = self.d_duration;
577        }
578        if self.v_s >= self.v_threshold {
579            self.v_s = self.v_reset;
580            1
581        } else {
582            0
583        }
584    }
585    pub fn reset(&mut self) {
586        self.v_s = -65.0;
587        self.v_d = -65.0;
588        self.d_active = false;
589        self.d_timer = 0.0;
590    }
591}
592impl Default for DendrifyNeuron {
593    fn default() -> Self {
594        Self::new()
595    }
596}
597
598/// Two-compartment LIF — soma + dendrite with history-dependent coupling.
599#[derive(Clone, Debug)]
600pub struct TwoCompartmentLIFNeuron {
601    pub v_s: f64,
602    pub v_d: f64,
603    pub v_rest: f64,
604    pub v_reset: f64,
605    pub theta: f64,
606    pub tau_s: f64,
607    pub tau_d: f64,
608    pub kappa: f64,
609    pub dt: f64,
610}
611
612impl TwoCompartmentLIFNeuron {
613    pub fn new() -> Self {
614        Self {
615            v_s: 0.0,
616            v_d: 0.0,
617            v_rest: 0.0,
618            v_reset: 0.0,
619            theta: 1.0,
620            tau_s: 2.0,
621            tau_d: 20.0,
622            kappa: 0.5,
623            dt: 1.0,
624        }
625    }
626    pub fn step(&mut self, i_soma: f64, i_dend: f64) -> i32 {
627        let alpha_s = (-self.dt / self.tau_s).exp();
628        let alpha_d = (-self.dt / self.tau_d).exp();
629        self.v_d = alpha_d * self.v_d + i_dend;
630        self.v_s = alpha_s * self.v_s + i_soma + self.kappa * self.v_d;
631        if self.v_s >= self.theta {
632            self.v_s = self.v_reset;
633            1
634        } else {
635            0
636        }
637    }
638    pub fn reset(&mut self) {
639        self.v_s = self.v_rest;
640        self.v_d = self.v_rest;
641    }
642}
643impl Default for TwoCompartmentLIFNeuron {
644    fn default() -> Self {
645        Self::new()
646    }
647}
648
649#[cfg(test)]
650mod tests {
651    use super::*;
652
653    #[test]
654    fn pr_fires() {
655        let mut n = PinskyRinzelNeuron::new();
656        let t: i32 = (0..5000).map(|_| n.step(5.0, 0.0)).sum();
657        assert!(t > 0);
658    }
659    #[test]
660    fn hay_fires() {
661        let mut n = HayL5PyramidalNeuron::new();
662        let t: i32 = (0..500).map(|_| n.step(20.0, 0.0)).sum();
663        assert!(t > 0);
664    }
665    #[test]
666    fn marder_fires() {
667        let mut n = MarderSTGNeuron::new();
668        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
669        assert!(t > 0);
670    }
671    #[test]
672    fn rall_fires() {
673        let mut n = RallCableNeuron::new(5);
674        let t: i32 = (0..500).map(|_| n.step(50.0)).sum();
675        assert!(t > 0);
676    }
677    #[test]
678    fn booth_fires() {
679        let mut n = BoothRinzelNeuron::new();
680        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
681        assert!(t > 0);
682    }
683    #[test]
684    fn dendrify_fires() {
685        let mut n = DendrifyNeuron::new();
686        let t: i32 = (0..2000).map(|_| n.step(50.0)).sum();
687        assert!(t > 0);
688    }
689    #[test]
690    fn tc_lif_fires() {
691        let mut n = TwoCompartmentLIFNeuron::new();
692        let t: i32 = (0..100).map(|_| n.step(0.5, 0.3)).sum();
693        assert!(t > 0);
694    }
695
696    // ── Multi-angle tests for multi-compartment models ──
697
698    // -- PinskyRinzel --
699    #[test]
700    fn pr_reset() {
701        let mut n = PinskyRinzelNeuron::new();
702        for _ in 0..100 {
703            n.step(5.0, 0.0);
704        }
705        n.reset();
706        assert!((n.v_s - (-60.0)).abs() < 1e-10);
707        assert!((n.v_d - (-60.0)).abs() < 1e-10);
708    }
709    #[test]
710    fn pr_bounded() {
711        let mut n = PinskyRinzelNeuron::new();
712        for _ in 0..5000 {
713            n.step(50.0, 0.0);
714        }
715        assert!(n.v_s.is_finite());
716    }
717    #[test]
718    fn pr_dendritic_input() {
719        let mut n = PinskyRinzelNeuron::new();
720        let _t: i32 = (0..5000).map(|_| n.step(0.0, 5.0)).sum();
721        // Dendritic input should also be able to drive spiking
722        assert!(n.v_d.is_finite());
723    }
724    #[test]
725    fn pr_nan_no_panic() {
726        PinskyRinzelNeuron::new().step(f64::NAN, 0.0);
727    }
728
729    // -- HayL5 --
730    #[test]
731    fn hay_reset() {
732        let mut n = HayL5PyramidalNeuron::new();
733        for _ in 0..100 {
734            n.step(20.0, 0.0);
735        }
736        n.reset();
737        assert!((n.v_s - (-75.0)).abs() < 1e-10);
738    }
739    #[test]
740    fn hay_bounded() {
741        let mut n = HayL5PyramidalNeuron::new();
742        for _ in 0..500 {
743            n.step(100.0, 0.0);
744        }
745        assert!(n.v_s.is_finite());
746    }
747    #[test]
748    fn hay_nan_no_panic() {
749        HayL5PyramidalNeuron::new().step(f64::NAN, 0.0);
750    }
751
752    // -- MarderSTG --
753    #[test]
754    fn marder_reset() {
755        let mut n = MarderSTGNeuron::new();
756        for _ in 0..100 {
757            n.step(5.0);
758        }
759        n.reset();
760        assert!((n.v - (-60.0)).abs() < 1e-10);
761    }
762    #[test]
763    fn marder_bounded() {
764        let mut n = MarderSTGNeuron::new();
765        for _ in 0..2000 {
766            n.step(50.0);
767        }
768        assert!(n.v.is_finite());
769    }
770    #[test]
771    fn marder_nan_no_panic() {
772        MarderSTGNeuron::new().step(f64::NAN);
773    }
774
775    // -- RallCable --
776    #[test]
777    fn rall_reset() {
778        let mut n = RallCableNeuron::new(5);
779        for _ in 0..100 {
780            n.step(50.0);
781        }
782        n.reset();
783        assert!(n.v.iter().all(|&x| (x - n.v_rest).abs() < 1e-10));
784    }
785    #[test]
786    fn rall_bounded() {
787        let mut n = RallCableNeuron::new(5);
788        for _ in 0..1000 {
789            n.step(500.0);
790        }
791        assert!(n.v.iter().all(|x| x.is_finite()));
792    }
793    #[test]
794    fn rall_nan_no_panic() {
795        RallCableNeuron::new(5).step(f64::NAN);
796    }
797
798    // -- BoothRinzel --
799    #[test]
800    fn booth_reset() {
801        let mut n = BoothRinzelNeuron::new();
802        for _ in 0..100 {
803            n.step(5.0);
804        }
805        n.reset();
806        assert!((n.vs - (-65.0)).abs() < 1e-10);
807    }
808    #[test]
809    fn booth_bounded() {
810        let mut n = BoothRinzelNeuron::new();
811        for _ in 0..2000 {
812            n.step(50.0);
813        }
814        assert!(n.vs.is_finite());
815    }
816    #[test]
817    fn booth_nan_no_panic() {
818        BoothRinzelNeuron::new().step(f64::NAN);
819    }
820
821    // -- Dendrify --
822    #[test]
823    fn dendrify_reset() {
824        let mut n = DendrifyNeuron::new();
825        for _ in 0..100 {
826            n.step(50.0);
827        }
828        n.reset();
829        assert!((n.v_s - (-65.0)).abs() < 1e-10);
830    }
831    #[test]
832    fn dendrify_bounded() {
833        let mut n = DendrifyNeuron::new();
834        for _ in 0..2000 {
835            n.step(200.0);
836        }
837        assert!(n.v_s.is_finite());
838    }
839    #[test]
840    fn dendrify_nan_no_panic() {
841        DendrifyNeuron::new().step(f64::NAN);
842    }
843
844    // -- TwoCompartmentLIF --
845    #[test]
846    fn tc_lif_reset() {
847        let mut n = TwoCompartmentLIFNeuron::new();
848        for _ in 0..50 {
849            n.step(0.5, 0.3);
850        }
851        n.reset();
852    }
853    #[test]
854    fn tc_lif_bounded() {
855        let mut n = TwoCompartmentLIFNeuron::new();
856        for _ in 0..1000 {
857            n.step(100.0, 100.0);
858        }
859        assert!(n.v_s.is_finite());
860    }
861    #[test]
862    fn tc_lif_nan_no_panic() {
863        TwoCompartmentLIFNeuron::new().step(f64::NAN, 0.0);
864    }
865}
866
867/// Dendritic NMDA spike model.
868///
869/// Captures the non-linear voltage-dependent Mg²⁺ block of NMDA receptors
870/// in dendritic branches. NMDA current has a sigmoidal voltage dependence:
871///
872///   I_NMDA = g_NMDA · B(V) · (V - E_NMDA)
873///   B(V) = 1 / (1 + [Mg²⁺]/3.57 · exp(-0.062 · V))
874///
875/// This enables coincidence detection: the dendrite only passes current
876/// when both presynaptic glutamate AND postsynaptic depolarisation are present.
877///
878/// Reference: Jahr & Stevens (1990), Schiller et al. (2000).
879#[derive(Clone, Debug)]
880pub struct DendriticNMDANeuron {
881    pub v_soma: f64,
882    pub v_dend: f64,
883    pub g_nmda: f64,
884    pub e_nmda: f64,
885    pub mg_conc: f64,
886    pub g_coupling: f64,
887    pub tau_soma: f64,
888    pub tau_dend: f64,
889    pub theta: f64,
890    pub dt: f64,
891}
892
893impl DendriticNMDANeuron {
894    pub fn new() -> Self {
895        Self {
896            v_soma: -65.0,
897            v_dend: -65.0,
898            g_nmda: 1.5,
899            e_nmda: 0.0,
900            mg_conc: 1.0,
901            g_coupling: 0.5,
902            tau_soma: 20.0,
903            tau_dend: 50.0,
904            theta: -50.0,
905            dt: 0.1,
906        }
907    }
908
909    /// Mg²⁺ block factor (Jahr & Stevens 1990).
910    fn mg_block(&self, v: f64) -> f64 {
911        1.0 / (1.0 + (self.mg_conc / 3.57) * (-0.062 * v).exp())
912    }
913
914    /// Step with somatic input and dendritic glutamate.
915    pub fn step(&mut self, i_soma: f64, glutamate: f64) -> i32 {
916        // Dendritic NMDA current with Mg²⁺ block.
917        let b = self.mg_block(self.v_dend);
918        let i_nmda = self.g_nmda * glutamate * b * (self.v_dend - self.e_nmda);
919
920        // Dendrite dynamics.
921        let dv_dend =
922            (-self.v_dend - 65.0 + i_nmda + self.g_coupling * (self.v_soma - self.v_dend))
923                / self.tau_dend;
924        self.v_dend += dv_dend * self.dt;
925
926        // Soma dynamics: receives dendritic current + direct input.
927        let i_dend_to_soma = self.g_coupling * (self.v_dend - self.v_soma);
928        let dv_soma = (-self.v_soma - 65.0 + i_soma + i_dend_to_soma) / self.tau_soma;
929        self.v_soma += dv_soma * self.dt;
930
931        if self.v_soma >= self.theta {
932            self.v_soma = -65.0;
933            1
934        } else {
935            0
936        }
937    }
938
939    pub fn reset(&mut self) {
940        self.v_soma = -65.0;
941        self.v_dend = -65.0;
942    }
943}
944
945impl Default for DendriticNMDANeuron {
946    fn default() -> Self {
947        Self::new()
948    }
949}
950
951/// Multi-compartment neuron (MCN) matching the Spiking-WM architecture.
952///
953/// Dual-dendrite model with basal and apical compartments. The apical dendrite
954/// gates how strongly basal information influences the soma, enabling
955/// nonlinear integration for long-term temporal memory in RL tasks.
956///
957/// Exact equations from arXiv:2503.00713 (Spiking-WM, PNAS 2025):
958///
959///   τ_b dV_b/dt = -V_b + x_b                                  (basal)
960///   τ_a dV_a/dt = -V_a + x_a                                  (apical)
961///   τ   dU/dt   = -U + σ(V_a)·[g_B/g_L·(V_b - U) + W_s·I]   (soma)
962///   S[t] = Θ(U[t] - V_th)                                     (spike)
963///   U[t] ← U[t]·(1 - S[t])                                    (soft reset)
964///
965/// Default parameters from Table II: τ = τ_a = τ_b = 2.0, g_B/g_L = 1.0,
966/// β = 1.0 (sigmoid steepness), V_th = 1.0.
967///
968/// Reference: Brain-Cog-Lab, arXiv:2503.00713, PNAS 2025.
969#[derive(Clone, Debug)]
970pub struct MulticompartmentMCNNeuron {
971    /// Somatic membrane potential.
972    pub u: f64,
973    /// Basal dendrite potential.
974    pub v_basal: f64,
975    /// Apical dendrite potential.
976    pub v_apical: f64,
977    /// Soma time constant.
978    pub tau: f64,
979    /// Basal dendrite time constant.
980    pub tau_b: f64,
981    /// Apical dendrite time constant.
982    pub tau_a: f64,
983    /// Basal-to-soma conductance ratio (g_B/g_L).
984    pub g_ratio: f64,
985    /// Sigmoid steepness for apical gating.
986    pub beta: f64,
987    /// Spike threshold.
988    pub v_th: f64,
989    /// Time step.
990    pub dt: f64,
991}
992
993impl MulticompartmentMCNNeuron {
994    pub fn new() -> Self {
995        Self {
996            u: 0.0,
997            v_basal: 0.0,
998            v_apical: 0.0,
999            tau: 2.0,
1000            tau_b: 2.0,
1001            tau_a: 2.0,
1002            g_ratio: 1.0,
1003            beta: 1.0,
1004            v_th: 1.0,
1005            dt: 1.0,
1006        }
1007    }
1008
1009    /// Sigmoid gating function σ(x) = 1/(1 + exp(-βx)).
1010    fn sigma(&self, x: f64) -> f64 {
1011        1.0 / (1.0 + (-self.beta * x).exp())
1012    }
1013
1014    /// Step with basal input (x_b), apical input (x_a), and direct somatic input.
1015    pub fn step_compartments(&mut self, x_basal: f64, x_apical: f64, i_soma: f64) -> i32 {
1016        // Basal dendrite: τ_b dV_b/dt = -V_b + x_b.
1017        let dv_b = (-self.v_basal + x_basal) / self.tau_b;
1018        self.v_basal += dv_b * self.dt;
1019
1020        // Apical dendrite: τ_a dV_a/dt = -V_a + x_a.
1021        let dv_a = (-self.v_apical + x_apical) / self.tau_a;
1022        self.v_apical += dv_a * self.dt;
1023
1024        // Soma: τ dU/dt = -U + σ(V_a)·[g_B/g_L·(V_b - U) + i_soma].
1025        let gate = self.sigma(self.v_apical);
1026        let du = (-self.u + gate * (self.g_ratio * (self.v_basal - self.u) + i_soma)) / self.tau;
1027        self.u += du * self.dt;
1028
1029        // Spike: S = Θ(U - V_th), reset: U ← U·(1 - S).
1030        if self.u >= self.v_th {
1031            self.u = 0.0;
1032            1
1033        } else {
1034            0
1035        }
1036    }
1037
1038    /// Simple step: input goes to basal dendrite only.
1039    pub fn step(&mut self, current: f64) -> i32 {
1040        self.step_compartments(current, 0.0, 0.0)
1041    }
1042
1043    pub fn reset(&mut self) {
1044        self.u = 0.0;
1045        self.v_basal = 0.0;
1046        self.v_apical = 0.0;
1047    }
1048}
1049
1050impl Default for MulticompartmentMCNNeuron {
1051    fn default() -> Self {
1052        Self::new()
1053    }
1054}
1055
1056/// Astrocyte-LIF hybrid unit with calcium wave feedback.
1057///
1058/// Models the tripartite synapse: a glial astrocyte monitors extracellular
1059/// glutamate from a paired LIF neuron and provides slow homeostatic feedback
1060/// via calcium-dependent gliotransmitter release.
1061///
1062///   dCa/dt = -Ca/τ_ca + δ · S_pre(t)        (calcium rise on presynaptic spike)
1063///   I_glio = g_glio · H(Ca - Ca_thresh)      (gliotransmitter release)
1064///   dV/dt = -(V - E_L)/τ_m + I_ext + I_glio  (LIF with glial feedback)
1065///
1066/// Reference: Perea, Navarrete & Araque, "Tripartite synapses" (2009).
1067#[derive(Clone, Debug)]
1068pub struct AstrocyteLIFNeuron {
1069    pub v: f64,
1070    pub ca: f64,
1071    pub tau_m: f64,
1072    pub tau_ca: f64,
1073    pub e_l: f64,
1074    pub theta: f64,
1075    pub v_reset: f64,
1076    pub ca_delta: f64,
1077    pub ca_thresh: f64,
1078    pub g_glio: f64,
1079    pub dt: f64,
1080}
1081
1082impl AstrocyteLIFNeuron {
1083    pub fn new() -> Self {
1084        Self {
1085            v: -65.0,
1086            ca: 0.0,
1087            tau_m: 20.0,
1088            tau_ca: 500.0,
1089            e_l: -65.0,
1090            theta: -50.0,
1091            v_reset: -65.0,
1092            ca_delta: 0.1,
1093            ca_thresh: 0.5,
1094            g_glio: 2.0,
1095            dt: 0.1,
1096        }
1097    }
1098
1099    /// Step with external current and presynaptic spike indicator.
1100    pub fn step_with_pre(&mut self, i_ext: f64, pre_spike: bool) -> i32 {
1101        // Astrocyte calcium dynamics.
1102        let dca = -self.ca / self.tau_ca
1103            + if pre_spike {
1104                self.ca_delta / self.dt
1105            } else {
1106                0.0
1107            };
1108        self.ca += dca * self.dt;
1109        self.ca = self.ca.max(0.0);
1110
1111        // Gliotransmitter release (Heaviside on calcium).
1112        let i_glio = if self.ca > self.ca_thresh {
1113            self.g_glio
1114        } else {
1115            0.0
1116        };
1117
1118        // LIF membrane dynamics with glial feedback.
1119        let dv = (-(self.v - self.e_l) + i_ext + i_glio) / self.tau_m;
1120        self.v += dv * self.dt;
1121
1122        if self.v >= self.theta {
1123            self.v = self.v_reset;
1124            1
1125        } else {
1126            0
1127        }
1128    }
1129
1130    /// Simple step (no presynaptic spike).
1131    pub fn step(&mut self, current: f64) -> i32 {
1132        self.step_with_pre(current, false)
1133    }
1134
1135    pub fn reset(&mut self) {
1136        self.v = self.e_l;
1137        self.ca = 0.0;
1138    }
1139}
1140
1141impl Default for AstrocyteLIFNeuron {
1142    fn default() -> Self {
1143        Self::new()
1144    }
1145}
1146
1147// ---- Tests for new multi-compartment / glial models ----
1148
1149#[cfg(test)]
1150mod gap_mc_tests {
1151    use super::*;
1152
1153    #[test]
1154    fn nmda_coincidence_detection() {
1155        let mut n = DendriticNMDANeuron::new();
1156        // Only soma input — dendrite contributes little.
1157        let mut spikes_soma_only = 0;
1158        for _ in 0..2000 {
1159            spikes_soma_only += n.step(8.0, 0.0);
1160        }
1161        n.reset();
1162        // Soma + glutamate — NMDA amplifies.
1163        let mut spikes_both = 0;
1164        for _ in 0..2000 {
1165            spikes_both += n.step(8.0, 1.0);
1166        }
1167        // Coincidence: both inputs together should fire more.
1168        assert!(
1169            spikes_both >= spikes_soma_only,
1170            "NMDA coincidence: both={spikes_both} must >= soma_only={spikes_soma_only}"
1171        );
1172    }
1173
1174    #[test]
1175    fn nmda_mg_block_voltage_dependent() {
1176        let n = DendriticNMDANeuron::new();
1177        let b_hyper = n.mg_block(-80.0);
1178        let b_depol = n.mg_block(-20.0);
1179        assert!(
1180            b_depol > b_hyper,
1181            "Mg block must relieve at depolarised potentials: B(-20)={b_depol:.3} > B(-80)={b_hyper:.3}"
1182        );
1183    }
1184
1185    #[test]
1186    fn nmda_zero_glutamate_no_nmda_current() {
1187        let mut n = DendriticNMDANeuron::new();
1188        let spikes: i32 = (0..500).map(|_| n.step(0.0, 0.0)).sum();
1189        assert_eq!(spikes, 0, "No input → no spikes");
1190    }
1191
1192    #[test]
1193    fn mcn_apical_gating() {
1194        // Without apical input, gate = σ(0) = 0.5, moderate drive.
1195        let mut n_no_apical = MulticompartmentMCNNeuron::new();
1196        let mut spikes_no = 0;
1197        for _ in 0..100 {
1198            spikes_no += n_no_apical.step_compartments(3.0, 0.0, 0.0);
1199        }
1200        // With strong apical input, gate ≈ 1.0, full basal→soma coupling.
1201        let mut n_apical = MulticompartmentMCNNeuron::new();
1202        let mut spikes_yes = 0;
1203        for _ in 0..100 {
1204            spikes_yes += n_apical.step_compartments(3.0, 5.0, 0.0);
1205        }
1206        assert!(
1207            spikes_yes >= spikes_no,
1208            "Apical gating should boost firing: apical={spikes_yes} >= none={spikes_no}"
1209        );
1210    }
1211
1212    #[test]
1213    fn mcn_basal_dendrite_memory() {
1214        // τ_b = 2.0, dt = 1.0: V_b decays by factor (1 - dt/τ) = 0.5 per step.
1215        let mut n = MulticompartmentMCNNeuron::new();
1216        n.step_compartments(5.0, 0.0, 0.0);
1217        let v_after = n.v_basal;
1218        n.step_compartments(0.0, 0.0, 0.0);
1219        let v_decay = n.v_basal;
1220        assert!(
1221            v_decay.abs() > 0.1 * v_after.abs(),
1222            "Basal dendrite retains memory: {v_decay:.3} vs {v_after:.3}"
1223        );
1224    }
1225
1226    #[test]
1227    fn mcn_reset_clears_all() {
1228        let mut n = MulticompartmentMCNNeuron::new();
1229        for _ in 0..50 {
1230            n.step(2.0);
1231        }
1232        n.reset();
1233        assert_eq!(n.u, 0.0);
1234        assert_eq!(n.v_basal, 0.0);
1235        assert_eq!(n.v_apical, 0.0);
1236    }
1237
1238    #[test]
1239    fn astrocyte_calcium_rises_on_pre_spikes() {
1240        let mut n = AstrocyteLIFNeuron::new();
1241        let ca_before = n.ca;
1242        for _ in 0..100 {
1243            n.step_with_pre(0.0, true);
1244        }
1245        assert!(
1246            n.ca > ca_before,
1247            "Calcium must rise with presynaptic spikes"
1248        );
1249    }
1250
1251    #[test]
1252    fn astrocyte_gliotransmitter_boosts_firing() {
1253        let mut n_no_glio = AstrocyteLIFNeuron::new();
1254        let mut n_glio = AstrocyteLIFNeuron::new();
1255
1256        let mut spikes_no = 0;
1257        let mut spikes_yes = 0;
1258        for _ in 0..5000 {
1259            spikes_no += n_no_glio.step_with_pre(10.0, false);
1260            spikes_yes += n_glio.step_with_pre(10.0, true); // pre spikes → Ca → glio
1261        }
1262        assert!(
1263            spikes_yes >= spikes_no,
1264            "Gliotransmitter should boost firing: with={spikes_yes} >= without={spikes_no}"
1265        );
1266    }
1267
1268    #[test]
1269    fn astrocyte_calcium_decays() {
1270        let mut n = AstrocyteLIFNeuron::new();
1271        // Build up calcium.
1272        for _ in 0..200 {
1273            n.step_with_pre(0.0, true);
1274        }
1275        let ca_peak = n.ca;
1276        // Let it decay.
1277        for _ in 0..5000 {
1278            n.step_with_pre(0.0, false);
1279        }
1280        assert!(
1281            n.ca < ca_peak * 0.5,
1282            "Calcium must decay: current={:.4} < peak={:.4}*0.5",
1283            n.ca,
1284            ca_peak
1285        );
1286    }
1287}
sc_neurocore_engine/neurons/multi_compartment.rs

sc_neurocore_engine/neurons/
multi_compartment.rs
