Skip to main content

sc_neurocore_engine/neurons/
multi_compartment.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
2// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
3// © Code 2020–2026 Miroslav Šotek. All rights reserved.
4// ORCID: 0009-0009-3560-0851
5// Contact: www.anulum.li | protoscience@anulum.li
6// SC-NeuroCore — Multi-compartment neuron models
7
8//! Multi-compartment neuron models.
9
10#[allow(unused_imports)]
11use super::biophysical::safe_rate;
12
13/// Pinsky-Rinzel 1994 — 2-compartment pyramidal cell.
14#[derive(Clone, Debug)]
15pub struct PinskyRinzelNeuron {
16    pub v_s: f64,
17    pub v_d: f64,
18    pub h: f64,
19    pub n: f64,
20    pub s: f64,
21    pub c: f64,
22    pub q: f64,
23    pub gc: f64,
24    pub p: f64,
25    pub g_na: f64,
26    pub g_kdr: f64,
27    pub g_ca: f64,
28    pub g_kahp: f64,
29    pub g_l: f64,
30    pub e_na: f64,
31    pub e_k: f64,
32    pub e_ca: f64,
33    pub e_l: f64,
34    pub dt: f64,
35    pub v_threshold: f64,
36}
37
38impl PinskyRinzelNeuron {
39    pub fn new() -> Self {
40        Self {
41            v_s: -60.0,
42            v_d: -60.0,
43            h: 0.9,
44            n: 0.1,
45            s: 0.0,
46            c: 0.0,
47            q: 0.0,
48            gc: 2.1,
49            p: 0.5,
50            g_na: 30.0,
51            g_kdr: 15.0,
52            g_ca: 10.0,
53            g_kahp: 0.8,
54            g_l: 0.1,
55            e_na: 60.0,
56            e_k: -75.0,
57            e_ca: 80.0,
58            e_l: -60.0,
59            dt: 0.02,
60            v_threshold: -20.0,
61        }
62    }
63    pub fn step(&mut self, current_soma: f64, current_dend: f64) -> i32 {
64        let v_prev = self.v_s;
65        let am = safe_rate(0.32, 54.0, self.v_s, 4.0, 8.0);
66        let bm = safe_rate(-0.28, 27.0, self.v_s, -5.0, 5.6);
67        let m_inf = am / (am + bm);
68        let ah = 0.128 * (-(self.v_s + 50.0) / 18.0).exp();
69        let bh = 4.0 / (1.0 + (-(self.v_s + 27.0) / 5.0).exp());
70        let an = safe_rate(0.032, 52.0, self.v_s, 5.0, 0.32);
71        let bn = 0.5 * (-(self.v_s + 57.0) / 40.0).exp();
72        let s_inf = 1.0 / (1.0 + (-(self.v_d + 20.0) / 9.0).exp());
73        let i_na = self.g_na * m_inf.powi(2) * self.h * (self.v_s - self.e_na);
74        let i_kdr = self.g_kdr * self.n.powi(2) * (self.v_s - self.e_k);
75        let i_ls = self.g_l * (self.v_s - self.e_l);
76        let i_ds = (self.gc / self.p) * (self.v_s - self.v_d);
77        let i_ca = self.g_ca * self.s.powi(2) * (self.v_d - self.e_ca);
78        let i_kahp = self.g_kahp * self.q * (self.v_d - self.e_k);
79        let i_ld = self.g_l * (self.v_d - self.e_l);
80        let i_sd = (self.gc / (1.0 - self.p)) * (self.v_d - self.v_s);
81        self.v_s += (-i_na - i_kdr - i_ls - i_ds + current_soma / self.p) * self.dt;
82        self.v_d += (-i_ca - i_kahp - i_ld - i_sd + current_dend / (1.0 - self.p)) * self.dt;
83        self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
84        self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
85        self.s += ((s_inf - self.s) / 5.0) * self.dt;
86        self.c = (self.c + (-0.13 * i_ca - 0.075 * self.c) * self.dt).max(0.0);
87        let q_inf = (self.c / (self.c + 2.0)).min(1.0);
88        self.q += ((q_inf - self.q) / 100.0) * self.dt;
89        if self.v_s >= self.v_threshold && v_prev < self.v_threshold {
90            1
91        } else {
92            0
93        }
94    }
95    pub fn reset(&mut self) {
96        self.v_s = -60.0;
97        self.v_d = -60.0;
98        self.h = 0.9;
99        self.n = 0.1;
100        self.s = 0.0;
101        self.c = 0.0;
102        self.q = 0.0;
103    }
104}
105impl Default for PinskyRinzelNeuron {
106    fn default() -> Self {
107        Self::new()
108    }
109}
110
111/// Hay et al. 2011 — Layer 5 thick-tufted pyramidal (3-compartment reduced).
112#[derive(Clone, Debug)]
113pub struct HayL5PyramidalNeuron {
114    pub v_s: f64,
115    pub v_t: f64,
116    pub v_a: f64,
117    pub h_na: f64,
118    pub n_k: f64,
119    pub m_ca: f64,
120    pub h_ca: f64,
121    pub m_ih: f64,
122    pub ca_a: f64,
123    pub g_na: f64,
124    pub g_k: f64,
125    pub g_l_s: f64,
126    pub g_ca_t: f64,
127    pub g_ih: f64,
128    pub g_l_t: f64,
129    pub g_ca_a: f64,
130    pub g_kca: f64,
131    pub g_l_a: f64,
132    pub g_st: f64,
133    pub g_ta: f64,
134    pub p_s: f64,
135    pub p_t: f64,
136    pub p_a: f64,
137    pub e_na: f64,
138    pub e_k: f64,
139    pub e_ca: f64,
140    pub e_ih: f64,
141    pub e_l: f64,
142    pub ca_decay: f64,
143    pub f_ca: f64,
144    pub c_m: f64,
145    pub dt: f64,
146    pub v_threshold: f64,
147}
148
149impl HayL5PyramidalNeuron {
150    pub fn new() -> Self {
151        Self {
152            v_s: -75.0,
153            v_t: -75.0,
154            v_a: -75.0,
155            h_na: 0.9,
156            n_k: 0.1,
157            m_ca: 0.0,
158            h_ca: 1.0,
159            m_ih: 0.0,
160            ca_a: 0.0001,
161            g_na: 300.0,
162            g_k: 40.0,
163            g_l_s: 0.03,
164            g_ca_t: 2.0,
165            g_ih: 0.02,
166            g_l_t: 0.03,
167            g_ca_a: 1.5,
168            g_kca: 2.5,
169            g_l_a: 0.03,
170            g_st: 1.5,
171            g_ta: 0.8,
172            p_s: 0.15,
173            p_t: 0.25,
174            p_a: 0.60,
175            e_na: 50.0,
176            e_k: -85.0,
177            e_ca: 140.0,
178            e_ih: -45.0,
179            e_l: -75.0,
180            ca_decay: 200.0,
181            f_ca: 0.0002,
182            c_m: 1.0,
183            dt: 0.025,
184            v_threshold: -30.0,
185        }
186    }
187    pub fn step(&mut self, current_soma: f64, current_tuft: f64) -> i32 {
188        let v_s_prev = self.v_s;
189        for _ in 0..4 {
190            let m_na_inf = 1.0 / (1.0 + (-(self.v_s + 38.0) / 7.0).exp());
191            let h_na_inf = 1.0 / (1.0 + ((self.v_s + 65.0) / 6.0).exp());
192            let n_k_inf = 1.0 / (1.0 + (-(self.v_s + 25.0) / 12.0).exp());
193            let tau_h = 0.5 + 14.0 / (1.0 + ((self.v_s + 35.0) / 10.0).exp());
194            let tau_n = 1.0 + 5.0 / (1.0 + ((self.v_s + 30.0) / 10.0).exp());
195            self.h_na += (h_na_inf - self.h_na) / tau_h * self.dt;
196            self.n_k += (n_k_inf - self.n_k) / tau_n * self.dt;
197            let i_na = self.g_na * m_na_inf.powi(3) * self.h_na * (self.v_s - self.e_na);
198            let i_k = self.g_k * self.n_k.powi(4) * (self.v_s - self.e_k);
199            let i_l_s = self.g_l_s * (self.v_s - self.e_l);
200            let i_st = self.g_st * (self.v_s - self.v_t) / self.p_s;
201            let m_ca_inf = 1.0 / (1.0 + (-(self.v_t + 27.0) / 7.0).exp());
202            let h_ca_inf = 1.0 / (1.0 + ((self.v_t + 52.0) / 5.0).exp());
203            let m_ih_inf = 1.0 / (1.0 + ((self.v_t + 75.0) / 5.5).exp());
204            self.m_ca += (m_ca_inf - self.m_ca) / 1.0 * self.dt;
205            self.h_ca += (h_ca_inf - self.h_ca) / 20.0 * self.dt;
206            self.m_ih += (m_ih_inf - self.m_ih) / 50.0 * self.dt;
207            let i_ca_t = self.g_ca_t * self.m_ca.powi(2) * self.h_ca * (self.v_t - self.e_ca);
208            let i_ih = self.g_ih * self.m_ih * (self.v_t - self.e_ih);
209            let i_l_t = self.g_l_t * (self.v_t - self.e_l);
210            let i_ts = self.g_st * (self.v_t - self.v_s) / self.p_t;
211            let i_ta = self.g_ta * (self.v_t - self.v_a) / self.p_t;
212            let m_ca_a_inf = 1.0 / (1.0 + (-(self.v_a + 30.0) / 5.0).exp());
213            let kca_act = self.ca_a / (self.ca_a + 0.001);
214            let i_ca_a = self.g_ca_a * m_ca_a_inf.powi(2) * (self.v_a - self.e_ca);
215            let i_kca = self.g_kca * kca_act * (self.v_a - self.e_k);
216            let i_l_a = self.g_l_a * (self.v_a - self.e_l);
217            let i_at = self.g_ta * (self.v_a - self.v_t) / self.p_a;
218            self.v_s += (-i_na - i_k - i_l_s - i_st + current_soma / self.p_s) / self.c_m * self.dt;
219            self.v_t += (-i_ca_t - i_ih - i_l_t - i_ts - i_ta) / self.c_m * self.dt;
220            self.v_a +=
221                (-i_ca_a - i_kca - i_l_a - i_at + current_tuft / self.p_a) / self.c_m * self.dt;
222            self.ca_a =
223                (self.ca_a + (-self.f_ca * i_ca_a - self.ca_a / self.ca_decay) * self.dt).max(0.0);
224        }
225        if self.v_s >= self.v_threshold && v_s_prev < self.v_threshold {
226            1
227        } else {
228            0
229        }
230    }
231    pub fn reset(&mut self) {
232        self.v_s = -75.0;
233        self.v_t = -75.0;
234        self.v_a = -75.0;
235        self.h_na = 0.9;
236        self.n_k = 0.1;
237        self.m_ca = 0.0;
238        self.h_ca = 1.0;
239        self.m_ih = 0.0;
240        self.ca_a = 0.0001;
241    }
242}
243impl Default for HayL5PyramidalNeuron {
244    fn default() -> Self {
245        Self::new()
246    }
247}
248
249/// Marder STG — stomatogastric ganglion LP neuron with 7 currents.
250#[derive(Clone, Debug)]
251pub struct MarderSTGNeuron {
252    pub v: f64,
253    pub m_na: f64,
254    pub h_na: f64,
255    pub m_cat: f64,
256    pub h_cat: f64,
257    pub m_cas: f64,
258    pub m_a: f64,
259    pub h_a: f64,
260    pub m_kd: f64,
261    pub m_h: f64,
262    pub ca: f64,
263    pub g_na: f64,
264    pub g_cat: f64,
265    pub g_cas: f64,
266    pub g_a: f64,
267    pub g_kd: f64,
268    pub g_kca: f64,
269    pub g_h: f64,
270    pub g_l: f64,
271    pub e_na: f64,
272    pub e_k: f64,
273    pub e_ca: f64,
274    pub e_h: f64,
275    pub e_l: f64,
276    pub dt: f64,
277    pub v_threshold: f64,
278}
279
280impl MarderSTGNeuron {
281    pub fn new() -> Self {
282        Self {
283            v: -60.0,
284            m_na: 0.0,
285            h_na: 0.9,
286            m_cat: 0.0,
287            h_cat: 0.9,
288            m_cas: 0.0,
289            m_a: 0.0,
290            h_a: 0.9,
291            m_kd: 0.0,
292            m_h: 0.0,
293            ca: 0.05,
294            g_na: 400.0,
295            g_cat: 2.5,
296            g_cas: 6.0,
297            g_a: 50.0,
298            g_kd: 100.0,
299            g_kca: 25.0,
300            g_h: 0.02,
301            g_l: 0.01,
302            e_na: 50.0,
303            e_k: -80.0,
304            e_ca: 120.0,
305            e_h: -20.0,
306            e_l: -50.0,
307            dt: 0.05,
308            v_threshold: -20.0,
309        }
310    }
311    pub fn step(&mut self, current: f64) -> i32 {
312        let v_prev = self.v;
313        let b = |v: f64, vh: f64, s: f64| 1.0 / (1.0 + (-(v - vh) / s).exp());
314        self.m_na += (b(self.v, -25.5, 5.29) - self.m_na) / 1.32 * self.dt;
315        self.h_na += (b(self.v, -48.9, -5.18) - self.h_na)
316            / (0.67 * (1.0 + ((self.v + 62.9) / -10.0).exp()) + 1.5)
317            * self.dt;
318        self.m_cat += (b(self.v, -27.1, 7.2) - self.m_cat) / 21.7 * self.dt;
319        self.h_cat += (b(self.v, -32.1, -5.5) - self.h_cat) / 105.0 * self.dt;
320        self.m_cas += (b(self.v, -33.0, 8.1) - self.m_cas) / 14.0 * self.dt;
321        self.m_a += (b(self.v, -27.2, 8.7) - self.m_a) / 11.6 * self.dt;
322        self.h_a += (b(self.v, -56.9, -4.9) - self.h_a) / 38.6 * self.dt;
323        self.m_kd += (b(self.v, -12.3, 11.8) - self.m_kd) / 7.2 * self.dt;
324        self.m_h += (b(self.v, -70.0, -6.0) - self.m_h) / 272.0 * self.dt;
325        let kca_act = (self.ca / (self.ca + 3.0)).min(1.0);
326        let i_na = self.g_na * self.m_na.powi(3) * self.h_na * (self.v - self.e_na);
327        let i_cat = self.g_cat * self.m_cat.powi(3) * self.h_cat * (self.v - self.e_ca);
328        let i_cas = self.g_cas * self.m_cas.powi(3) * (self.v - self.e_ca);
329        let i_a = self.g_a * self.m_a.powi(3) * self.h_a * (self.v - self.e_k);
330        let i_kd = self.g_kd * self.m_kd.powi(4) * (self.v - self.e_k);
331        let i_kca = self.g_kca * kca_act.powi(4) * (self.v - self.e_k);
332        let i_h = self.g_h * self.m_h * (self.v - self.e_h);
333        let i_l = self.g_l * (self.v - self.e_l);
334        self.v += (-i_na - i_cat - i_cas - i_a - i_kd - i_kca - i_h - i_l + current) * self.dt;
335        let i_ca_total = i_cat + i_cas;
336        self.ca = (self.ca + (-0.0001 * i_ca_total - 0.01 * self.ca) * self.dt).max(0.0);
337        if self.v >= self.v_threshold && v_prev < self.v_threshold {
338            1
339        } else {
340            0
341        }
342    }
343    pub fn reset(&mut self) {
344        self.v = -60.0;
345        self.m_na = 0.0;
346        self.h_na = 0.9;
347        self.m_cat = 0.0;
348        self.h_cat = 0.9;
349        self.m_cas = 0.0;
350        self.m_a = 0.0;
351        self.h_a = 0.9;
352        self.m_kd = 0.0;
353        self.m_h = 0.0;
354        self.ca = 0.05;
355    }
356}
357impl Default for MarderSTGNeuron {
358    fn default() -> Self {
359        Self::new()
360    }
361}
362
363/// Rall cable — N-compartment passive dendrite model. Rall 1964.
364#[derive(Clone, Debug)]
365pub struct RallCableNeuron {
366    pub v: Vec<f64>,
367    pub n_comp: usize,
368    pub tau_m: f64,
369    pub v_rest: f64,
370    pub g_ratio: f64,
371    pub v_threshold: f64,
372    pub v_reset: f64,
373    pub dt: f64,
374}
375
376impl RallCableNeuron {
377    pub fn new(n_comp: usize) -> Self {
378        Self {
379            v: vec![-65.0; n_comp],
380            n_comp,
381            tau_m: 20.0,
382            v_rest: -65.0,
383            g_ratio: 0.5,
384            v_threshold: -50.0,
385            v_reset: -65.0,
386            dt: 0.1,
387        }
388    }
389    pub fn step(&mut self, current: f64) -> i32 {
390        let n = self.n_comp;
391        let mut dv = vec![0.0; n];
392        for i in 0..n {
393            let leak = -(self.v[i] - self.v_rest) / self.tau_m;
394            let left = if i > 0 {
395                self.g_ratio * (self.v[i - 1] - self.v[i])
396            } else {
397                0.0
398            };
399            let right = if i + 1 < n {
400                self.g_ratio * (self.v[i + 1] - self.v[i])
401            } else {
402                0.0
403            };
404            let inj = if i == n - 1 { current } else { 0.0 };
405            dv[i] = (leak + left + right + inj) * self.dt;
406        }
407        for i in 0..n {
408            self.v[i] += dv[i];
409        }
410        if self.v[0] >= self.v_threshold {
411            self.v[0] = self.v_reset;
412            1
413        } else {
414            0
415        }
416    }
417    pub fn reset(&mut self) {
418        self.v.fill(self.v_rest);
419    }
420}
421
422/// Booth-Rinzel — 2-compartment motoneuron with bistability. Booth et al. 1997.
423#[derive(Clone, Debug)]
424pub struct BoothRinzelNeuron {
425    pub vs: f64,
426    pub vd: f64,
427    pub h: f64,
428    pub n: f64,
429    pub q: f64,
430    pub ca: f64,
431    pub p: f64,
432    pub gc: f64,
433    pub g_na: f64,
434    pub g_k: f64,
435    pub g_ca: f64,
436    pub g_kca: f64,
437    pub g_l: f64,
438    pub e_na: f64,
439    pub e_k: f64,
440    pub e_ca: f64,
441    pub e_l: f64,
442    pub dt: f64,
443    pub v_threshold: f64,
444}
445
446impl BoothRinzelNeuron {
447    pub fn new() -> Self {
448        Self {
449            vs: -65.0,
450            vd: -65.0,
451            h: 0.9,
452            n: 0.0,
453            q: 0.0,
454            ca: 0.0,
455            p: 0.5,
456            gc: 0.1,
457            g_na: 120.0,
458            g_k: 100.0,
459            g_ca: 14.0,
460            g_kca: 5.0,
461            g_l: 0.51,
462            e_na: 55.0,
463            e_k: -80.0,
464            e_ca: 80.0,
465            e_l: -60.0,
466            dt: 0.025,
467            v_threshold: -20.0,
468        }
469    }
470    pub fn step(&mut self, current: f64) -> i32 {
471        let v_prev = self.vs;
472        for _ in 0..4 {
473            let m_inf = 1.0 / (1.0 + (-(self.vs + 30.0) / 9.0).exp());
474            let h_inf = 1.0 / (1.0 + ((self.vs + 45.0) / 7.0).exp());
475            let n_inf = 1.0 / (1.0 + (-(self.vs + 30.0) / 10.0).exp());
476            let m_ca_inf = 1.0 / (1.0 + (-(self.vd + 20.0) / 9.0).exp());
477            let q_inf = (self.ca / (self.ca + 2.0)).min(1.0);
478            self.h += (h_inf - self.h) / 1.0 * self.dt;
479            self.n += (n_inf - self.n) / 3.0 * self.dt;
480            self.q += (q_inf - self.q) / 100.0 * self.dt;
481            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.vs - self.e_na);
482            let i_k = self.g_k * self.n.powi(4) * (self.vs - self.e_k);
483            let i_ls = self.g_l * (self.vs - self.e_l);
484            let i_sd = (self.gc / self.p) * (self.vs - self.vd);
485            let i_ca = self.g_ca * m_ca_inf.powi(2) * (self.vd - self.e_ca);
486            let i_kca = self.g_kca * self.q * (self.vd - self.e_k);
487            let i_ld = self.g_l * (self.vd - self.e_l);
488            let i_ds = (self.gc / (1.0 - self.p)) * (self.vd - self.vs);
489            self.vs += (-i_na - i_k - i_ls - i_sd + current / self.p) * self.dt;
490            self.vd += (-i_ca - i_kca - i_ld - i_ds) * self.dt;
491            self.ca = (self.ca + (-0.13 * i_ca - 0.075 * self.ca) * self.dt).max(0.0);
492        }
493        if self.vs >= self.v_threshold && v_prev < self.v_threshold {
494            1
495        } else {
496            0
497        }
498    }
499    pub fn reset(&mut self) {
500        self.vs = -65.0;
501        self.vd = -65.0;
502        self.h = 0.9;
503        self.n = 0.0;
504        self.q = 0.0;
505        self.ca = 0.0;
506    }
507}
508impl Default for BoothRinzelNeuron {
509    fn default() -> Self {
510        Self::new()
511    }
512}
513
514/// Dendrify — two-compartment with active dendritic spike (NMDA-like plateau).
515#[derive(Clone, Debug)]
516pub struct DendrifyNeuron {
517    pub v_s: f64,
518    pub v_d: f64,
519    pub d_active: bool,
520    pub d_timer: f64,
521    pub tau_s: f64,
522    pub tau_d: f64,
523    pub g_c: f64,
524    pub d_threshold: f64,
525    pub d_amplitude: f64,
526    pub d_duration: f64,
527    pub v_rest: f64,
528    pub v_threshold: f64,
529    pub v_reset: f64,
530    pub dt: f64,
531}
532
533impl DendrifyNeuron {
534    pub fn new() -> Self {
535        Self {
536            v_s: -65.0,
537            v_d: -65.0,
538            d_active: false,
539            d_timer: 0.0,
540            tau_s: 10.0,
541            tau_d: 20.0,
542            g_c: 0.8,
543            d_threshold: -35.0,
544            d_amplitude: 30.0,
545            d_duration: 10.0,
546            v_rest: -65.0,
547            v_threshold: -50.0,
548            v_reset: -65.0,
549            dt: 0.1,
550        }
551    }
552    pub fn step(&mut self, current: f64) -> i32 {
553        let d_input = if self.d_active { self.d_amplitude } else { 0.0 };
554        self.v_d += (-(self.v_d - self.v_rest) + current + d_input
555            - self.g_c * (self.v_d - self.v_s))
556            / self.tau_d
557            * self.dt;
558        self.v_s +=
559            (-(self.v_s - self.v_rest) + self.g_c * (self.v_d - self.v_s)) / self.tau_s * self.dt;
560        if self.d_active {
561            self.d_timer -= self.dt;
562            if self.d_timer <= 0.0 {
563                self.d_active = false;
564            }
565        } else if self.v_d >= self.d_threshold {
566            self.d_active = true;
567            self.d_timer = self.d_duration;
568        }
569        if self.v_s >= self.v_threshold {
570            self.v_s = self.v_reset;
571            1
572        } else {
573            0
574        }
575    }
576    pub fn reset(&mut self) {
577        self.v_s = -65.0;
578        self.v_d = -65.0;
579        self.d_active = false;
580        self.d_timer = 0.0;
581    }
582}
583impl Default for DendrifyNeuron {
584    fn default() -> Self {
585        Self::new()
586    }
587}
588
589/// Two-compartment LIF — soma + dendrite with history-dependent coupling.
590#[derive(Clone, Debug)]
591pub struct TwoCompartmentLIFNeuron {
592    pub v_s: f64,
593    pub v_d: f64,
594    pub v_rest: f64,
595    pub v_reset: f64,
596    pub theta: f64,
597    pub tau_s: f64,
598    pub tau_d: f64,
599    pub kappa: f64,
600    pub dt: f64,
601}
602
603impl TwoCompartmentLIFNeuron {
604    pub fn new() -> Self {
605        Self {
606            v_s: 0.0,
607            v_d: 0.0,
608            v_rest: 0.0,
609            v_reset: 0.0,
610            theta: 1.0,
611            tau_s: 2.0,
612            tau_d: 20.0,
613            kappa: 0.5,
614            dt: 1.0,
615        }
616    }
617    pub fn step(&mut self, i_soma: f64, i_dend: f64) -> i32 {
618        let alpha_s = (-self.dt / self.tau_s).exp();
619        let alpha_d = (-self.dt / self.tau_d).exp();
620        self.v_d = alpha_d * self.v_d + i_dend;
621        self.v_s = alpha_s * self.v_s + i_soma + self.kappa * self.v_d;
622        if self.v_s >= self.theta {
623            self.v_s = self.v_reset;
624            1
625        } else {
626            0
627        }
628    }
629    pub fn reset(&mut self) {
630        self.v_s = self.v_rest;
631        self.v_d = self.v_rest;
632    }
633}
634impl Default for TwoCompartmentLIFNeuron {
635    fn default() -> Self {
636        Self::new()
637    }
638}
639
640#[cfg(test)]
641mod tests {
642    use super::*;
643
644    #[test]
645    fn pr_fires() {
646        let mut n = PinskyRinzelNeuron::new();
647        let t: i32 = (0..5000).map(|_| n.step(5.0, 0.0)).sum();
648        assert!(t > 0);
649    }
650    #[test]
651    fn hay_fires() {
652        let mut n = HayL5PyramidalNeuron::new();
653        let t: i32 = (0..500).map(|_| n.step(20.0, 0.0)).sum();
654        assert!(t > 0);
655    }
656    #[test]
657    fn marder_fires() {
658        let mut n = MarderSTGNeuron::new();
659        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
660        assert!(t > 0);
661    }
662    #[test]
663    fn rall_fires() {
664        let mut n = RallCableNeuron::new(5);
665        let t: i32 = (0..500).map(|_| n.step(50.0)).sum();
666        assert!(t > 0);
667    }
668    #[test]
669    fn booth_fires() {
670        let mut n = BoothRinzelNeuron::new();
671        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
672        assert!(t > 0);
673    }
674    #[test]
675    fn dendrify_fires() {
676        let mut n = DendrifyNeuron::new();
677        let t: i32 = (0..2000).map(|_| n.step(50.0)).sum();
678        assert!(t > 0);
679    }
680    #[test]
681    fn tc_lif_fires() {
682        let mut n = TwoCompartmentLIFNeuron::new();
683        let t: i32 = (0..100).map(|_| n.step(0.5, 0.3)).sum();
684        assert!(t > 0);
685    }
686}