Skip to main content

sc_neurocore_engine/neurons/
simple_spiking.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 — Two-dimensional and higher spiking neuron models
7
8//! Two-dimensional and higher spiking neuron models.
9
10use rand::{RngExt, SeedableRng};
11use rand_xoshiro::Xoshiro256PlusPlus;
12
13/// FitzHugh-Nagumo 1961 — 2D qualitative spike model.
14#[derive(Clone, Debug)]
15pub struct FitzHughNagumoNeuron {
16    pub v: f64,
17    pub w: f64,
18    pub a: f64,
19    pub b: f64,
20    pub epsilon: f64,
21    pub dt: f64,
22    pub v_threshold: f64,
23}
24
25impl FitzHughNagumoNeuron {
26    pub fn new() -> Self {
27        Self {
28            v: -1.0,
29            w: -0.5,
30            a: 0.7,
31            b: 0.8,
32            epsilon: 0.08,
33            dt: 0.1,
34            v_threshold: 1.0,
35        }
36    }
37    pub fn step(&mut self, current: f64) -> i32 {
38        let v_prev = self.v;
39        self.v += (self.v - self.v.powi(3) / 3.0 - self.w + current) * self.dt;
40        self.w += self.epsilon * (self.v + self.a - self.b * self.w) * self.dt;
41        if self.v >= self.v_threshold && v_prev < self.v_threshold {
42            1
43        } else {
44            0
45        }
46    }
47    pub fn reset(&mut self) {
48        self.v = -1.0;
49        self.w = -0.5;
50    }
51}
52impl Default for FitzHughNagumoNeuron {
53    fn default() -> Self {
54        Self::new()
55    }
56}
57
58/// Morris-Lecar 1981 — barnacle muscle fiber.
59#[derive(Clone, Debug)]
60pub struct MorrisLecarNeuron {
61    pub v: f64,
62    pub w: f64,
63    pub c_m: f64,
64    pub g_ca: f64,
65    pub g_k: f64,
66    pub g_l: f64,
67    pub e_ca: f64,
68    pub e_k: f64,
69    pub e_l: f64,
70    pub v1: f64,
71    pub v2: f64,
72    pub v3: f64,
73    pub v4: f64,
74    pub phi: f64,
75    pub dt: f64,
76    pub v_threshold: f64,
77}
78
79impl MorrisLecarNeuron {
80    pub fn new() -> Self {
81        Self {
82            v: -60.0,
83            w: 0.0,
84            c_m: 20.0,
85            g_ca: 4.0,
86            g_k: 8.0,
87            g_l: 2.0,
88            e_ca: 120.0,
89            e_k: -84.0,
90            e_l: -60.0,
91            v1: -1.2,
92            v2: 18.0,
93            v3: 12.0,
94            v4: 17.4,
95            phi: 1.0 / 15.0,
96            dt: 0.1,
97            v_threshold: 0.0,
98        }
99    }
100    pub fn step(&mut self, current: f64) -> i32 {
101        let v_prev = self.v;
102        let m_inf = 0.5 * (1.0 + ((self.v - self.v1) / self.v2).tanh());
103        let w_inf = 0.5 * (1.0 + ((self.v - self.v3) / self.v4).tanh());
104        let lam = self.phi * ((self.v - self.v3) / (2.0 * self.v4)).cosh();
105        let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
106        let i_k = self.g_k * self.w * (self.v - self.e_k);
107        let i_l = self.g_l * (self.v - self.e_l);
108        self.v += (-i_ca - i_k - i_l + current) / self.c_m * self.dt;
109        self.w += lam * (w_inf - self.w) * self.dt;
110        if self.v >= self.v_threshold && v_prev < self.v_threshold {
111            1
112        } else {
113            0
114        }
115    }
116    pub fn reset(&mut self) {
117        self.v = -60.0;
118        self.w = 0.0;
119    }
120}
121impl Default for MorrisLecarNeuron {
122    fn default() -> Self {
123        Self::new()
124    }
125}
126
127/// Hindmarsh-Rose 1984 — 3D chaotic bursting model.
128#[derive(Clone, Debug)]
129pub struct HindmarshRoseNeuron {
130    pub x: f64,
131    pub y: f64,
132    pub z: f64,
133    pub b: f64,
134    pub r: f64,
135    pub s: f64,
136    pub x_rest: f64,
137    pub dt: f64,
138    pub x_threshold: f64,
139}
140
141impl HindmarshRoseNeuron {
142    pub fn new() -> Self {
143        Self {
144            x: -1.6,
145            y: -10.0,
146            z: 2.0,
147            b: 3.0,
148            r: 0.001,
149            s: 4.0,
150            x_rest: -1.6,
151            dt: 0.1,
152            x_threshold: 1.0,
153        }
154    }
155    pub fn step(&mut self, current: f64) -> i32 {
156        let x_prev = self.x;
157        let dx = (self.y - self.x.powi(3) + self.b * self.x.powi(2) - self.z + current) * self.dt;
158        let dy = (1.0 - 5.0 * self.x.powi(2) - self.y) * self.dt;
159        let dz = self.r * (self.s * (self.x - self.x_rest) - self.z) * self.dt;
160        self.x += dx;
161        self.y += dy;
162        self.z += dz;
163        if self.x >= self.x_threshold && x_prev < self.x_threshold {
164            1
165        } else {
166            0
167        }
168    }
169    pub fn reset(&mut self) {
170        self.x = -1.6;
171        self.y = -10.0;
172        self.z = 2.0;
173    }
174}
175impl Default for HindmarshRoseNeuron {
176    fn default() -> Self {
177        Self::new()
178    }
179}
180
181/// Resonate-and-Fire — damped oscillator with threshold. Izhikevich 2001.
182#[derive(Clone, Debug)]
183pub struct ResonateAndFireNeuron {
184    pub x: f64,
185    pub y: f64,
186    pub b: f64,
187    pub omega: f64,
188    pub threshold: f64,
189    pub dt: f64,
190}
191
192impl ResonateAndFireNeuron {
193    pub fn new() -> Self {
194        Self {
195            x: 0.0,
196            y: 0.0,
197            b: -0.1,
198            omega: 1.0,
199            threshold: 1.0,
200            dt: 0.05,
201        }
202    }
203    pub fn step(&mut self, current: f64) -> i32 {
204        self.x += (self.b * self.x - self.omega * self.y + current) * self.dt;
205        self.y += (self.omega * self.x + self.b * self.y) * self.dt;
206        let r = (self.x * self.x + self.y * self.y).sqrt();
207        if r >= self.threshold {
208            self.x = 0.0;
209            self.y = 0.0;
210            1
211        } else {
212            0
213        }
214    }
215    pub fn reset(&mut self) {
216        self.x = 0.0;
217        self.y = 0.0;
218    }
219}
220impl Default for ResonateAndFireNeuron {
221    fn default() -> Self {
222        Self::new()
223    }
224}
225
226/// FitzHugh-Rinzel — 3D extension with slow bursting variable.
227#[derive(Clone, Debug)]
228pub struct FitzHughRinzelNeuron {
229    pub v: f64,
230    pub w: f64,
231    pub y: f64,
232    pub a: f64,
233    pub b: f64,
234    pub c: f64,
235    pub d: f64,
236    pub delta: f64,
237    pub mu: f64,
238    pub dt: f64,
239    pub v_threshold: f64,
240}
241
242impl FitzHughRinzelNeuron {
243    pub fn new() -> Self {
244        Self {
245            v: -1.0,
246            w: -0.5,
247            y: 0.0,
248            a: 0.7,
249            b: 0.8,
250            c: -0.775,
251            d: 1.0,
252            delta: 0.08,
253            mu: 0.0001,
254            dt: 0.1,
255            v_threshold: 1.0,
256        }
257    }
258    pub fn step(&mut self, current: f64) -> i32 {
259        let v_prev = self.v;
260        self.v += (self.v - self.v.powi(3) / 3.0 - self.w + self.y + current) * self.dt;
261        self.w += self.delta * (self.a + self.v - self.b * self.w) * self.dt;
262        self.y += self.mu * (self.c - self.v - self.d * self.y) * self.dt;
263        if self.v >= self.v_threshold && v_prev < self.v_threshold {
264            1
265        } else {
266            0
267        }
268    }
269    pub fn reset(&mut self) {
270        self.v = -1.0;
271        self.w = -0.5;
272        self.y = 0.0;
273    }
274}
275impl Default for FitzHughRinzelNeuron {
276    fn default() -> Self {
277        Self::new()
278    }
279}
280
281/// McKean 1970 — piecewise-linear FHN variant.
282#[derive(Clone, Debug)]
283pub struct McKeanNeuron {
284    pub v: f64,
285    pub w: f64,
286    pub a: f64,
287    pub epsilon: f64,
288    pub gamma: f64,
289    pub dt: f64,
290    pub v_peak: f64,
291}
292
293impl McKeanNeuron {
294    pub fn new() -> Self {
295        Self {
296            v: 0.0,
297            w: 0.0,
298            a: 0.25,
299            epsilon: 0.01,
300            gamma: 0.5,
301            dt: 0.1,
302            v_peak: 0.8,
303        }
304    }
305    fn f_v(&self, v: f64) -> f64 {
306        let half_a = self.a / 2.0;
307        let mid = (1.0 + self.a) / 2.0;
308        if v < half_a {
309            -v
310        } else if v < mid {
311            v - self.a
312        } else {
313            1.0 - v
314        }
315    }
316    pub fn step(&mut self, current: f64) -> i32 {
317        let v_prev = self.v;
318        self.v += (self.f_v(self.v) - self.w + current) * self.dt;
319        self.w += self.epsilon * (self.v - self.gamma * self.w) * self.dt;
320        if self.v >= self.v_peak && v_prev < self.v_peak {
321            1
322        } else {
323            0
324        }
325    }
326    pub fn reset(&mut self) {
327        self.v = 0.0;
328        self.w = 0.0;
329    }
330}
331impl Default for McKeanNeuron {
332    fn default() -> Self {
333        Self::new()
334    }
335}
336
337/// Terman-Wang oscillator — segmentation by oscillatory correlation.
338#[derive(Clone, Debug)]
339pub struct TermanWangOscillator {
340    pub v: f64,
341    pub w: f64,
342    pub alpha: f64,
343    pub beta: f64,
344    pub epsilon: f64,
345    pub rho: f64,
346    pub dt: f64,
347    pub v_peak: f64,
348}
349
350impl TermanWangOscillator {
351    pub fn new() -> Self {
352        Self {
353            v: -1.5,
354            w: -0.5,
355            alpha: 3.0,
356            beta: 0.2,
357            epsilon: 0.02,
358            rho: 0.0,
359            dt: 0.05,
360            v_peak: 1.5,
361        }
362    }
363    pub fn step(&mut self, current: f64) -> i32 {
364        let v_prev = self.v;
365        let f = 3.0 * self.v - self.v.powi(3) + 2.0;
366        let g = self.alpha * (1.0 + (self.v / self.beta).tanh());
367        self.v += (f - self.w + current + self.rho) * self.dt;
368        self.w += self.epsilon * (g - self.w) * self.dt;
369        if self.v >= self.v_peak && v_prev < self.v_peak {
370            1
371        } else {
372            0
373        }
374    }
375    pub fn reset(&mut self) {
376        self.v = -1.5;
377        self.w = -0.5;
378    }
379}
380impl Default for TermanWangOscillator {
381    fn default() -> Self {
382        Self::new()
383    }
384}
385
386/// Benda-Herz 2003 — firing-rate adaptation via subtractive feedback.
387#[derive(Clone, Debug)]
388pub struct BendaHerzNeuron {
389    pub a: f64,
390    pub f_max: f64,
391    pub beta: f64,
392    pub i_half: f64,
393    pub tau_a: f64,
394    pub delta_a: f64,
395    pub dt: f64,
396    rng: Xoshiro256PlusPlus,
397}
398
399impl BendaHerzNeuron {
400    pub fn new(seed: u64) -> Self {
401        Self {
402            a: 0.0,
403            f_max: 200.0,
404            beta: 0.1,
405            i_half: 5.0,
406            tau_a: 100.0,
407            delta_a: 0.5,
408            dt: 1.0,
409            rng: Xoshiro256PlusPlus::seed_from_u64(seed),
410        }
411    }
412    pub fn step(&mut self, current: f64) -> i32 {
413        let x = current - self.a;
414        let rate = self.f_max / (1.0 + (-self.beta * (x - self.i_half)).exp());
415        self.a += (-self.a / self.tau_a + self.delta_a * rate) * self.dt;
416        let p = rate * self.dt / 1000.0;
417        if self.rng.random::<f64>() < p {
418            1
419        } else {
420            0
421        }
422    }
423    pub fn reset(&mut self) {
424        self.a = 0.0;
425    }
426}
427
428/// Alpha-synapse LIF — separate excitatory/inhibitory exponential synapses.
429#[derive(Clone, Debug)]
430pub struct AlphaNeuron {
431    pub v: f64,
432    pub i_exc: f64,
433    pub i_inh: f64,
434    pub v_rest: f64,
435    pub v_threshold: f64,
436    pub tau_v: f64,
437    pub tau_exc: f64,
438    pub tau_inh: f64,
439    pub dt: f64,
440}
441
442impl AlphaNeuron {
443    pub fn new() -> Self {
444        Self {
445            v: 0.0,
446            i_exc: 0.0,
447            i_inh: 0.0,
448            v_rest: 0.0,
449            v_threshold: 1.0,
450            tau_v: 20.0,
451            tau_exc: 5.0,
452            tau_inh: 10.0,
453            dt: 1.0,
454        }
455    }
456    pub fn step(&mut self, exc_current: f64, inh_current: f64) -> i32 {
457        self.i_exc += (-self.i_exc / self.tau_exc + exc_current) * self.dt;
458        self.i_inh += (-self.i_inh / self.tau_inh + inh_current) * self.dt;
459        self.v += (-(self.v - self.v_rest) + self.i_exc - self.i_inh) / self.tau_v * self.dt;
460        if self.v >= self.v_threshold {
461            self.v = self.v_rest;
462            1
463        } else {
464            0
465        }
466    }
467    pub fn reset(&mut self) {
468        self.v = self.v_rest;
469        self.i_exc = 0.0;
470        self.i_inh = 0.0;
471    }
472}
473impl Default for AlphaNeuron {
474    fn default() -> Self {
475        Self::new()
476    }
477}
478
479/// Conductance-based LIF (COBA). Brette et al. 2007.
480#[derive(Clone, Debug)]
481pub struct COBALIFNeuron {
482    pub v: f64,
483    pub g_e: f64,
484    pub g_i: f64,
485    pub c_m: f64,
486    pub g_l: f64,
487    pub e_l: f64,
488    pub e_e: f64,
489    pub e_i: f64,
490    pub tau_e: f64,
491    pub tau_i: f64,
492    pub v_threshold: f64,
493    pub v_reset: f64,
494    pub dt: f64,
495}
496
497impl COBALIFNeuron {
498    pub fn new() -> Self {
499        Self {
500            v: -65.0,
501            g_e: 0.0,
502            g_i: 0.0,
503            c_m: 200.0,
504            g_l: 10.0,
505            e_l: -65.0,
506            e_e: 0.0,
507            e_i: -80.0,
508            tau_e: 5.0,
509            tau_i: 10.0,
510            v_threshold: -50.0,
511            v_reset: -65.0,
512            dt: 0.1,
513        }
514    }
515    pub fn step(&mut self, current: f64, delta_ge: f64, delta_gi: f64) -> i32 {
516        self.g_e += delta_ge;
517        self.g_i += delta_gi;
518        let i_syn = self.g_e * (self.v - self.e_e) + self.g_i * (self.v - self.e_i);
519        self.v += (-self.g_l * (self.v - self.e_l) - i_syn + current) / self.c_m * self.dt;
520        self.g_e *= (-self.dt / self.tau_e).exp();
521        self.g_i *= (-self.dt / self.tau_i).exp();
522        if self.v >= self.v_threshold {
523            self.v = self.v_reset;
524            1
525        } else {
526            0
527        }
528    }
529    pub fn reset(&mut self) {
530        self.v = self.e_l;
531        self.g_e = 0.0;
532        self.g_i = 0.0;
533    }
534}
535impl Default for COBALIFNeuron {
536    fn default() -> Self {
537        Self::new()
538    }
539}
540
541/// Gutkin-Ermentrout — reduced cortical neuron with Type-I excitability.
542#[derive(Clone, Debug)]
543pub struct GutkinErmentroutNeuron {
544    pub v: f64,
545    pub n: f64,
546    pub g_na: f64,
547    pub g_k: f64,
548    pub g_l: f64,
549    pub e_na: f64,
550    pub e_k: f64,
551    pub e_l: f64,
552    pub dt: f64,
553    pub v_threshold: f64,
554}
555
556impl GutkinErmentroutNeuron {
557    pub fn new() -> Self {
558        Self {
559            v: -65.0,
560            n: 0.1,
561            g_na: 20.0,
562            g_k: 10.0,
563            g_l: 8.0,
564            e_na: 60.0,
565            e_k: -90.0,
566            e_l: -80.0,
567            dt: 0.05,
568            v_threshold: -20.0,
569        }
570    }
571    pub fn step(&mut self, current: f64) -> i32 {
572        let v_prev = self.v;
573        let m_inf = 1.0 / (1.0 + (-(self.v + 20.0) / 15.0).exp());
574        let n_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 5.0).exp());
575        self.n += (n_inf - self.n) * self.dt;
576        let i_na = self.g_na * m_inf * (self.v - self.e_na);
577        let i_k = self.g_k * self.n * (self.v - self.e_k);
578        let i_l = self.g_l * (self.v - self.e_l);
579        self.v += (-i_na - i_k - i_l + current) * self.dt;
580        if self.v >= self.v_threshold && v_prev < self.v_threshold {
581            1
582        } else {
583            0
584        }
585    }
586    pub fn reset(&mut self) {
587        self.v = -65.0;
588        self.n = 0.1;
589    }
590}
591impl Default for GutkinErmentroutNeuron {
592    fn default() -> Self {
593        Self::new()
594    }
595}
596
597/// Wilson HR — simplified cortical neuron. Wilson 1999.
598#[derive(Clone, Debug)]
599pub struct WilsonHRNeuron {
600    pub v: f64,
601    pub r: f64,
602    pub tau_r: f64,
603    pub v_peak: f64,
604    pub dt: f64,
605}
606
607impl WilsonHRNeuron {
608    pub fn new() -> Self {
609        Self {
610            v: -0.7,
611            r: 0.1,
612            tau_r: 1.9,
613            v_peak: 0.4,
614            dt: 0.05,
615        }
616    }
617    pub fn step(&mut self, current: f64) -> i32 {
618        let v_prev = self.v;
619        let poly = -(17.81 + 47.71 * self.v + 32.63 * self.v * self.v) * (self.v - 0.55);
620        let syn = -26.0 * self.r * (self.v + 0.92);
621        self.v += (poly + syn + current) * self.dt;
622        self.r += (-self.r + 1.35 * self.v + 1.03) / self.tau_r * self.dt;
623        if self.v >= self.v_peak && v_prev < self.v_peak {
624            1
625        } else {
626            0
627        }
628    }
629    pub fn reset(&mut self) {
630        self.v = -0.7;
631        self.r = 0.1;
632    }
633}
634impl Default for WilsonHRNeuron {
635    fn default() -> Self {
636        Self::new()
637    }
638}
639
640/// Chay 1985 — pancreatic beta cell with Ca-dependent K.
641#[derive(Clone, Debug)]
642pub struct ChayNeuron {
643    pub v: f64,
644    pub n: f64,
645    pub ca: f64,
646    pub g_ca: f64,
647    pub g_k: f64,
648    pub g_kca: f64,
649    pub g_l: f64,
650    pub e_ca: f64,
651    pub e_k: f64,
652    pub e_l: f64,
653    pub rho: f64,
654    pub alpha_ca: f64,
655    pub k_ca: f64,
656    pub dt: f64,
657    pub v_threshold: f64,
658}
659
660impl ChayNeuron {
661    pub fn new() -> Self {
662        Self {
663            v: -50.0,
664            n: 0.1,
665            ca: 0.1,
666            g_ca: 25.0,
667            g_k: 1400.0,
668            g_kca: 12.0,
669            g_l: 7.0,
670            e_ca: 100.0,
671            e_k: -75.0,
672            e_l: -40.0,
673            rho: 0.00015,
674            alpha_ca: 0.002,
675            k_ca: 0.04,
676            dt: 0.02,
677            v_threshold: -20.0,
678        }
679    }
680    pub fn step(&mut self, current: f64) -> i32 {
681        let v_prev = self.v;
682        let m_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 8.0).exp());
683        let n_inf = 1.0 / (1.0 + (-(self.v + 18.0) / 14.0).exp());
684        let d = (self.v + 18.0).abs().max(0.01);
685        let tau_n = 1.0 / (0.01 * d);
686        let kca_act = self.ca / (self.ca + 1.0);
687        let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
688        let i_k = self.g_k * self.n * (self.v - self.e_k);
689        let i_kca = self.g_kca * kca_act * (self.v - self.e_k);
690        let i_l = self.g_l * (self.v - self.e_l);
691        self.v += (-i_ca - i_k - i_kca - i_l + current) * self.dt;
692        self.n += (n_inf - self.n) / tau_n.max(0.01) * self.dt;
693        self.ca =
694            (self.ca + self.rho * (-self.alpha_ca * i_ca - self.k_ca * self.ca) * self.dt).max(0.0);
695        if self.v >= self.v_threshold && v_prev < self.v_threshold {
696            1
697        } else {
698            0
699        }
700    }
701    pub fn reset(&mut self) {
702        self.v = -50.0;
703        self.n = 0.1;
704        self.ca = 0.1;
705    }
706}
707impl Default for ChayNeuron {
708    fn default() -> Self {
709        Self::new()
710    }
711}
712
713/// Chay-Keizer — modified beta cell with inactivating Ca current.
714#[derive(Clone, Debug)]
715pub struct ChayKeizerNeuron {
716    pub v: f64,
717    pub n: f64,
718    pub ca: f64,
719    pub g_ca: f64,
720    pub g_k: f64,
721    pub g_kca: f64,
722    pub g_l: f64,
723    pub e_ca: f64,
724    pub e_k: f64,
725    pub e_l: f64,
726    pub k_d: f64,
727    pub f_ca: f64,
728    pub k_ca: f64,
729    pub dt: f64,
730    pub v_threshold: f64,
731}
732
733impl ChayKeizerNeuron {
734    pub fn new() -> Self {
735        Self {
736            v: -50.0,
737            n: 0.01,
738            ca: 0.1,
739            g_ca: 20.0,
740            g_k: 25.0,
741            g_kca: 12.0,
742            g_l: 0.1,
743            e_ca: 100.0,
744            e_k: -75.0,
745            e_l: -40.0,
746            k_d: 1.0,
747            f_ca: 0.004,
748            k_ca: 0.03,
749            dt: 0.02,
750            v_threshold: -20.0,
751        }
752    }
753    pub fn step(&mut self, current: f64) -> i32 {
754        let v_prev = self.v;
755        let m_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 8.0).exp());
756        let n_inf = 1.0 / (1.0 + (-(self.v + 18.0) / 7.0).exp());
757        let tau_n = 20.0;
758        let q_kca = self.ca / (self.ca + self.k_d);
759        let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
760        let i_k = self.g_k * self.n * (self.v - self.e_k);
761        let i_kca = self.g_kca * q_kca * (self.v - self.e_k);
762        let i_l = self.g_l * (self.v - self.e_l);
763        self.v += (-i_ca - i_k - i_kca - i_l + current) * self.dt;
764        self.n += (n_inf - self.n) / tau_n * self.dt;
765        self.ca = (self.ca + (-self.f_ca * i_ca - self.k_ca * self.ca) * self.dt).max(0.0);
766        if self.v >= self.v_threshold && v_prev < self.v_threshold {
767            1
768        } else {
769            0
770        }
771    }
772    pub fn reset(&mut self) {
773        self.v = -50.0;
774        self.n = 0.01;
775        self.ca = 0.1;
776    }
777}
778impl Default for ChayKeizerNeuron {
779    fn default() -> Self {
780        Self::new()
781    }
782}
783
784/// Sherman-Rinzel-Keizer 1988 — pancreatic beta cell (reduced).
785#[derive(Clone, Debug)]
786pub struct ShermanRinzelKeizerNeuron {
787    pub v: f64,
788    pub n: f64,
789    pub s: f64,
790    pub g_ca: f64,
791    pub g_k: f64,
792    pub g_s: f64,
793    pub e_ca: f64,
794    pub e_k: f64,
795    pub tau_s: f64,
796    pub dt: f64,
797    pub v_threshold: f64,
798}
799
800impl ShermanRinzelKeizerNeuron {
801    pub fn new() -> Self {
802        Self {
803            v: -50.0,
804            n: 0.1,
805            s: 0.1,
806            g_ca: 3.6,
807            g_k: 10.0,
808            g_s: 4.0,
809            e_ca: 25.0,
810            e_k: -75.0,
811            tau_s: 5000.0,
812            dt: 0.5,
813            v_threshold: -20.0,
814        }
815    }
816    pub fn step(&mut self, current: f64) -> i32 {
817        let v_prev = self.v;
818        let m_inf = 1.0 / (1.0 + (-(self.v + 20.0) / 12.0).exp());
819        let n_inf = 1.0 / (1.0 + (-(self.v + 16.0) / 5.6).exp());
820        let s_inf = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
821        let tau_n = 9.09;
822        let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
823        let i_k = self.g_k * self.n * (self.v - self.e_k);
824        let i_s = self.g_s * self.s * (self.v - self.e_k);
825        self.v += (-i_ca - i_k - i_s + current) * self.dt / 5.3;
826        self.n += (n_inf - self.n) / tau_n * self.dt;
827        self.s += (s_inf - self.s) / self.tau_s * self.dt;
828        if self.v >= self.v_threshold && v_prev < self.v_threshold {
829            1
830        } else {
831            0
832        }
833    }
834    pub fn reset(&mut self) {
835        self.v = -50.0;
836        self.n = 0.1;
837        self.s = 0.1;
838    }
839}
840impl Default for ShermanRinzelKeizerNeuron {
841    fn default() -> Self {
842        Self::new()
843    }
844}
845
846/// Butera pre-Bötzinger respiratory neuron. Butera et al. 1999.
847#[derive(Clone, Debug)]
848pub struct ButeraRespiratoryNeuron {
849    pub v: f64,
850    pub n: f64,
851    pub h_nap: f64,
852    pub g_na: f64,
853    pub g_nap: f64,
854    pub g_k: f64,
855    pub g_l: f64,
856    pub e_na: f64,
857    pub e_k: f64,
858    pub e_l: f64,
859    pub tau_h: f64,
860    pub dt: f64,
861    pub v_threshold: f64,
862}
863
864impl ButeraRespiratoryNeuron {
865    pub fn new() -> Self {
866        Self {
867            v: -50.0,
868            n: 0.01,
869            h_nap: 0.5,
870            g_na: 28.0,
871            g_nap: 2.8,
872            g_k: 11.2,
873            g_l: 2.8,
874            e_na: 50.0,
875            e_k: -85.0,
876            e_l: -65.0,
877            tau_h: 10000.0,
878            dt: 0.1,
879            v_threshold: -20.0,
880        }
881    }
882    pub fn step(&mut self, current: f64) -> i32 {
883        let v_prev = self.v;
884        let m_na = 1.0 / (1.0 + (-(self.v + 34.0) / 5.0).exp());
885        let n_inf = 1.0 / (1.0 + (-(self.v + 29.0) / 4.0).exp());
886        let m_nap = 1.0 / (1.0 + (-(self.v + 40.0) / 6.0).exp());
887        let h_nap_inf = 1.0 / (1.0 + ((self.v + 48.0) / 6.0).exp());
888        let tau_n = 10.0;
889        let h_na = 1.0 / (1.0 + ((self.v + 48.0) / 4.0).exp());
890        let i_na = self.g_na * m_na.powi(3) * h_na * (self.v - self.e_na);
891        let i_nap = self.g_nap * m_nap * self.h_nap * (self.v - self.e_na);
892        let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
893        let i_l = self.g_l * (self.v - self.e_l);
894        self.v += (-i_na - i_nap - i_k - i_l + current) * self.dt;
895        self.n += (n_inf - self.n) / tau_n * self.dt;
896        self.h_nap += (h_nap_inf - self.h_nap) / self.tau_h * self.dt;
897        if self.v >= self.v_threshold && v_prev < self.v_threshold {
898            1
899        } else {
900            0
901        }
902    }
903    pub fn reset(&mut self) {
904        self.v = -50.0;
905        self.n = 0.01;
906        self.h_nap = 0.5;
907    }
908}
909impl Default for ButeraRespiratoryNeuron {
910    fn default() -> Self {
911        Self::new()
912    }
913}
914
915/// e-prop ALIF — adaptive LIF for eligibility-propagation learning. Bellec et al. 2020.
916#[derive(Clone, Debug)]
917pub struct EPropALIFNeuron {
918    pub v: f64,
919    pub a: f64,
920    pub e_trace: f64,
921    pub alpha_m: f64,
922    pub alpha_a: f64,
923    pub v_threshold_base: f64,
924    pub beta: f64,
925}
926
927impl EPropALIFNeuron {
928    pub fn new(tau_m: f64, tau_a: f64, dt: f64) -> Self {
929        Self {
930            v: 0.0,
931            a: 0.0,
932            e_trace: 0.0,
933            alpha_m: (-dt / tau_m).exp(),
934            alpha_a: (-dt / tau_a).exp(),
935            v_threshold_base: 1.0,
936            beta: 0.07,
937        }
938    }
939    pub fn step(&mut self, current: f64) -> i32 {
940        self.v = self.alpha_m * self.v + current;
941        let threshold = self.v_threshold_base + self.beta * self.a;
942        let psi = ((1.0 - (self.v - threshold).abs()) * 0.3).max(0.0);
943        self.e_trace = self.alpha_a * self.e_trace + psi;
944        if self.v >= threshold {
945            self.v = 0.0;
946            self.a = self.alpha_a * self.a + 1.0;
947            1
948        } else {
949            self.a *= self.alpha_a;
950            0
951        }
952    }
953    pub fn reset(&mut self) {
954        self.v = 0.0;
955        self.a = 0.0;
956        self.e_trace = 0.0;
957    }
958}
959
960impl Default for EPropALIFNeuron {
961    fn default() -> Self {
962        Self::new(20.0, 200.0, 1.0)
963    }
964}
965
966/// SuperSpike — LIF with surrogate gradient trace. Zenke & Ganguli 2018.
967#[derive(Clone, Debug)]
968pub struct SuperSpikeNeuron {
969    pub v: f64,
970    pub trace: f64,
971    pub alpha_m: f64,
972    pub alpha_e: f64,
973    pub v_threshold: f64,
974    pub v_reset: f64,
975    pub beta_sg: f64,
976}
977
978impl SuperSpikeNeuron {
979    pub fn new(tau_m: f64, tau_e: f64, dt: f64) -> Self {
980        Self {
981            v: 0.0,
982            trace: 0.0,
983            alpha_m: (-dt / tau_m).exp(),
984            alpha_e: (-dt / tau_e).exp(),
985            v_threshold: 1.0,
986            v_reset: 0.0,
987            beta_sg: 10.0,
988        }
989    }
990    pub fn step(&mut self, current: f64) -> i32 {
991        self.v = self.alpha_m * self.v + current;
992        let sg = 1.0 / (self.beta_sg * (self.v - self.v_threshold).abs() + 1.0).powi(2);
993        self.trace = self.alpha_e * self.trace + sg;
994        if self.v >= self.v_threshold {
995            self.v = self.v_reset;
996            1
997        } else {
998            0
999        }
1000    }
1001    pub fn reset(&mut self) {
1002        self.v = 0.0;
1003        self.trace = 0.0;
1004    }
1005}
1006
1007impl Default for SuperSpikeNeuron {
1008    fn default() -> Self {
1009        Self::new(10.0, 10.0, 1.0)
1010    }
1011}
1012
1013/// Learnable Neuron Model (LNM) — parameterised activation + decay.
1014#[derive(Clone, Debug)]
1015pub struct LearnableNeuronModel {
1016    pub v: f64,
1017    pub alpha: f64,
1018    pub beta: f64,
1019    pub gamma: f64,
1020    pub v_threshold: f64,
1021    pub f_slope: f64,
1022    pub f_shift: f64,
1023}
1024
1025impl LearnableNeuronModel {
1026    pub fn new() -> Self {
1027        Self {
1028            v: 0.0,
1029            alpha: 0.9,
1030            beta: 0.1,
1031            gamma: 0.05,
1032            v_threshold: 1.0,
1033            f_slope: 5.0,
1034            f_shift: 0.5,
1035        }
1036    }
1037    pub fn step(&mut self, current: f64) -> i32 {
1038        let f_v = 1.0 / (1.0 + (-(self.f_slope * (self.v - self.f_shift))).exp());
1039        self.v = self.alpha * self.v + self.beta * current + self.gamma * f_v;
1040        if self.v >= self.v_threshold {
1041            self.v = 0.0;
1042            1
1043        } else {
1044            0
1045        }
1046    }
1047    pub fn reset(&mut self) {
1048        self.v = 0.0;
1049    }
1050}
1051impl Default for LearnableNeuronModel {
1052    fn default() -> Self {
1053        Self::new()
1054    }
1055}
1056
1057/// Pernarowski 1994 — simplified beta cell burster (3 ODE).
1058#[derive(Clone, Debug)]
1059pub struct PernarowskiNeuron {
1060    pub v: f64,
1061    pub w: f64,
1062    pub z: f64,
1063    pub alpha: f64,
1064    pub beta: f64,
1065    pub eps1: f64,
1066    pub eps2: f64,
1067    pub gamma: f64,
1068    pub dt: f64,
1069    pub v_threshold: f64,
1070}
1071
1072impl PernarowskiNeuron {
1073    pub fn new() -> Self {
1074        Self {
1075            v: -1.0,
1076            w: 0.0,
1077            z: 0.0,
1078            alpha: 0.1,
1079            beta: 0.5,
1080            eps1: 0.1,
1081            eps2: 0.001,
1082            gamma: 0.5,
1083            dt: 0.1,
1084            v_threshold: 0.5,
1085        }
1086    }
1087    pub fn step(&mut self, current: f64) -> i32 {
1088        let v_prev = self.v;
1089        let f_v = self.v - self.v.powi(3) / 3.0;
1090        self.v += (f_v - self.w - self.z + current) * self.dt;
1091        self.w += self.eps1 * (self.v - self.gamma * self.w + self.alpha) * self.dt;
1092        self.z += self.eps2 * (self.beta * (self.v + 0.7) - self.z) * self.dt;
1093        if self.v >= self.v_threshold && v_prev < self.v_threshold {
1094            1
1095        } else {
1096            0
1097        }
1098    }
1099    pub fn reset(&mut self) {
1100        self.v = -1.0;
1101        self.w = 0.0;
1102        self.z = 0.0;
1103    }
1104}
1105impl Default for PernarowskiNeuron {
1106    fn default() -> Self {
1107        Self::new()
1108    }
1109}
1110
1111#[cfg(test)]
1112mod tests {
1113    use super::*;
1114
1115    #[test]
1116    fn fhn_fires() {
1117        let mut n = FitzHughNagumoNeuron::new();
1118        let t: i32 = (0..2000).map(|_| n.step(1.0)).sum();
1119        assert!(t > 0);
1120    }
1121    #[test]
1122    fn morris_lecar_fires() {
1123        let mut n = MorrisLecarNeuron::new();
1124        let t: i32 = (0..2000).map(|_| n.step(100.0)).sum();
1125        assert!(t > 0);
1126    }
1127    #[test]
1128    fn hr_fires() {
1129        let mut n = HindmarshRoseNeuron::new();
1130        let t: i32 = (0..2000).map(|_| n.step(3.0)).sum();
1131        assert!(t > 0);
1132    }
1133    #[test]
1134    fn rnf_fires() {
1135        let mut n = ResonateAndFireNeuron::new();
1136        let t: i32 = (0..5000).map(|_| n.step(3.0)).sum();
1137        assert!(t > 0);
1138    }
1139    #[test]
1140    fn fhr_fires() {
1141        let mut n = FitzHughRinzelNeuron::new();
1142        let t: i32 = (0..2000).map(|_| n.step(1.0)).sum();
1143        assert!(t > 0);
1144    }
1145    #[test]
1146    fn mckean_fires() {
1147        let mut n = McKeanNeuron::new();
1148        let t: i32 = (0..2000).map(|_| n.step(0.5)).sum();
1149        assert!(t > 0);
1150    }
1151    #[test]
1152    fn tw_fires() {
1153        let mut n = TermanWangOscillator::new();
1154        let t: i32 = (0..2000).map(|_| n.step(0.5)).sum();
1155        assert!(t > 0);
1156    }
1157    #[test]
1158    fn benda_herz_fires() {
1159        let mut n = BendaHerzNeuron::new(42);
1160        let t: i32 = (0..10000).map(|_| n.step(20.0)).sum();
1161        assert!(t > 0);
1162    }
1163    #[test]
1164    fn alpha_fires() {
1165        let mut n = AlphaNeuron::new();
1166        let t: i32 = (0..100).map(|_| n.step(0.5, 0.0)).sum();
1167        assert!(t > 0);
1168    }
1169    #[test]
1170    fn coba_fires() {
1171        let mut n = COBALIFNeuron::new();
1172        let t: i32 = (0..2000).map(|_| n.step(500.0, 0.0, 0.0)).sum();
1173        assert!(t > 0);
1174    }
1175    #[test]
1176    fn gutkin_fires() {
1177        let mut n = GutkinErmentroutNeuron::new();
1178        let t: i32 = (0..2000).map(|_| n.step(15.0)).sum();
1179        assert!(t > 0);
1180    }
1181    #[test]
1182    fn wilson_hr_fires() {
1183        let mut n = WilsonHRNeuron::new();
1184        let t: i32 = (0..2000).map(|_| n.step(0.5)).sum();
1185        assert!(t > 0);
1186    }
1187    #[test]
1188    fn chay_fires() {
1189        let mut n = ChayNeuron::new();
1190        let t: i32 = (0..5000).map(|_| n.step(20.0)).sum();
1191        assert!(t > 0);
1192    }
1193    #[test]
1194    fn chay_keizer_fires() {
1195        let mut n = ChayKeizerNeuron::new();
1196        let t: i32 = (0..5000).map(|_| n.step(10.0)).sum();
1197        assert!(t > 0);
1198    }
1199    #[test]
1200    fn srk_fires() {
1201        let mut n = ShermanRinzelKeizerNeuron::new();
1202        let t: i32 = (0..5000).map(|_| n.step(5.0)).sum();
1203        assert!(t > 0);
1204    }
1205    #[test]
1206    fn butera_fires() {
1207        let mut n = ButeraRespiratoryNeuron::new();
1208        let t: i32 = (0..20000).map(|_| n.step(50.0)).sum();
1209        assert!(t > 0);
1210    }
1211    #[test]
1212    fn eprop_fires() {
1213        let mut n = EPropALIFNeuron::default();
1214        let t: i32 = (0..50).map(|_| n.step(0.5)).sum();
1215        assert!(t > 0);
1216    }
1217    #[test]
1218    fn superspike_fires() {
1219        let mut n = SuperSpikeNeuron::default();
1220        let t: i32 = (0..50).map(|_| n.step(0.5)).sum();
1221        assert!(t > 0);
1222    }
1223    #[test]
1224    fn lnm_fires() {
1225        let mut n = LearnableNeuronModel::new();
1226        let t: i32 = (0..50).map(|_| n.step(2.0)).sum();
1227        assert!(t > 0);
1228    }
1229    #[test]
1230    fn pernarowski_fires() {
1231        let mut n = PernarowskiNeuron::new();
1232        let t: i32 = (0..2000).map(|_| n.step(1.0)).sum();
1233        assert!(t > 0);
1234    }
1235}