sc_neurocore_engine/neurons/
biophysical.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 — Hodgkin-Huxley-type conductance-based neuron models
7
8//! Hodgkin-Huxley-type conductance-based neuron models.
9
10use rand::{RngExt, SeedableRng};
11use rand_xoshiro::Xoshiro256PlusPlus;
12
13// ── Helper: safe alpha/beta kinetics avoiding division-by-zero ──
14
15pub fn safe_rate(
16    number_factor: f64,
17    v_offset: f64,
18    v: f64,
19    denom_scale: f64,
20    fallback: f64,
21) -> f64 {
22    let d = v + v_offset;
23    if d.abs() < 1e-7 {
24        fallback
25    } else {
26        number_factor * d / (1.0 - (-d / denom_scale).exp())
27    }
28}
29
30/// Hodgkin-Huxley 1952 — 4-ODE ion channel model.
31#[derive(Clone, Debug)]
32pub struct HodgkinHuxleyNeuron {
33    pub v: f64,
34    pub m: f64,
35    pub h: f64,
36    pub n: f64,
37    pub c_m: f64,
38    pub g_na: f64,
39    pub g_k: f64,
40    pub g_l: f64,
41    pub e_na: f64,
42    pub e_k: f64,
43    pub e_l: f64,
44    pub dt: f64,
45    pub v_threshold: f64,
46}
47
48impl HodgkinHuxleyNeuron {
49    pub fn new() -> Self {
50        Self {
51            v: -65.0,
52            m: 0.05,
53            h: 0.6,
54            n: 0.32,
55            c_m: 1.0,
56            g_na: 120.0,
57            g_k: 36.0,
58            g_l: 0.3,
59            e_na: 50.0,
60            e_k: -77.0,
61            e_l: -54.4,
62            dt: 0.01,
63            v_threshold: 0.0,
64        }
65    }
66    pub fn step(&mut self, current: f64) -> i32 {
67        let v_prev = self.v;
68        let steps = (1.0 / self.dt) as usize;
69        for _ in 0..steps {
70            let am = safe_rate(0.1, 40.0, self.v, 10.0, 1.0);
71            let bm = 4.0 * (-(self.v + 65.0) / 18.0).exp();
72            let ah = 0.07 * (-(self.v + 65.0) / 20.0).exp();
73            let bh = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
74            let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
75            let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
76            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
77            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
78            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
79            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
80            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
81            let i_l = self.g_l * (self.v - self.e_l);
82            self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
83        }
84        if self.v >= self.v_threshold && v_prev < self.v_threshold {
85            1
86        } else {
87            0
88        }
89    }
90    pub fn reset(&mut self) {
91        self.v = -65.0;
92        self.m = 0.05;
93        self.h = 0.6;
94        self.n = 0.32;
95    }
96}
97impl Default for HodgkinHuxleyNeuron {
98    fn default() -> Self {
99        Self::new()
100    }
101}
102
103/// Traub-Miles — hippocampal pyramidal HH variant. Traub et al. 1991.
104#[derive(Clone, Debug)]
105pub struct TraubMilesNeuron {
106    pub v: f64,
107    pub m: f64,
108    pub h: f64,
109    pub n: f64,
110    pub g_na: f64,
111    pub g_k: f64,
112    pub g_l: f64,
113    pub e_na: f64,
114    pub e_k: f64,
115    pub e_l: f64,
116    pub dt: f64,
117    pub v_threshold: f64,
118}
119
120impl TraubMilesNeuron {
121    pub fn new() -> Self {
122        Self {
123            v: -67.0,
124            m: 0.05,
125            h: 0.6,
126            n: 0.3,
127            g_na: 100.0,
128            g_k: 80.0,
129            g_l: 0.1,
130            e_na: 50.0,
131            e_k: -100.0,
132            e_l: -67.0,
133            dt: 0.01,
134            v_threshold: -20.0,
135        }
136    }
137    pub fn step(&mut self, current: f64) -> i32 {
138        let v_prev = self.v;
139        for _ in 0..10 {
140            let am = safe_rate(0.32, 54.0, self.v, 4.0, 8.0);
141            let bm = safe_rate(-0.28, 27.0, self.v, -5.0, 5.6);
142            let ah = 0.128 * (-(self.v + 50.0) / 18.0).exp();
143            let bh = 4.0 / (1.0 + (-(self.v + 27.0) / 5.0).exp());
144            let an = safe_rate(0.032, 52.0, self.v, 5.0, 0.32);
145            let bn = 0.5 * (-(self.v + 57.0) / 40.0).exp();
146            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
147            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
148            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
149            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
150            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
151            let i_l = self.g_l * (self.v - self.e_l);
152            self.v += (-i_na - i_k - i_l + current) * self.dt;
153        }
154        if self.v >= self.v_threshold && v_prev < self.v_threshold {
155            1
156        } else {
157            0
158        }
159    }
160    pub fn reset(&mut self) {
161        self.v = -67.0;
162        self.m = 0.05;
163        self.h = 0.6;
164        self.n = 0.3;
165    }
166}
167impl Default for TraubMilesNeuron {
168    fn default() -> Self {
169        Self::new()
170    }
171}
172
173/// Wang-Buzsaki — fast-spiking GABAergic interneuron. Wang & Buzsáki 1996.
174#[derive(Clone, Debug)]
175pub struct WangBuzsakiNeuron {
176    pub v: f64,
177    pub h: f64,
178    pub n: f64,
179    pub g_na: f64,
180    pub g_k: f64,
181    pub g_l: f64,
182    pub e_na: f64,
183    pub e_k: f64,
184    pub e_l: f64,
185    pub c_m: f64,
186    pub phi: f64,
187    pub dt: f64,
188    pub v_threshold: f64,
189}
190
191impl WangBuzsakiNeuron {
192    pub fn new() -> Self {
193        Self {
194            v: -65.0,
195            h: 0.8,
196            n: 0.1,
197            g_na: 35.0,
198            g_k: 9.0,
199            g_l: 0.1,
200            e_na: 55.0,
201            e_k: -90.0,
202            e_l: -65.0,
203            c_m: 1.0,
204            phi: 5.0,
205            dt: 0.01,
206            v_threshold: -20.0,
207        }
208    }
209    pub fn step(&mut self, current: f64) -> i32 {
210        let v_prev = self.v;
211        for _ in 0..10 {
212            let am = safe_rate(0.1, 35.0, self.v, 10.0, 1.0);
213            let bm = 4.0 * (-(self.v + 60.0) / 18.0).exp();
214            let m_inf = am / (am + bm);
215            let ah = 0.07 * (-(self.v + 58.0) / 20.0).exp();
216            let bh = 1.0 / (1.0 + (-(self.v + 28.0) / 10.0).exp());
217            let an = safe_rate(0.01, 34.0, self.v, 10.0, 0.1);
218            let bn = 0.125 * (-(self.v + 44.0) / 80.0).exp();
219            self.h += self.phi * (ah * (1.0 - self.h) - bh * self.h) * self.dt;
220            self.n += self.phi * (an * (1.0 - self.n) - bn * self.n) * self.dt;
221            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.v - self.e_na);
222            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
223            let i_l = self.g_l * (self.v - self.e_l);
224            self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
225        }
226        if self.v >= self.v_threshold && v_prev < self.v_threshold {
227            1
228        } else {
229            0
230        }
231    }
232    pub fn reset(&mut self) {
233        self.v = -65.0;
234        self.h = 0.8;
235        self.n = 0.1;
236    }
237}
238impl Default for WangBuzsakiNeuron {
239    fn default() -> Self {
240        Self::new()
241    }
242}
243
244/// Connor-Stevens — A-type K current for delay tuning. Connor et al. 1977.
245#[derive(Clone, Debug)]
246pub struct ConnorStevensNeuron {
247    pub v: f64,
248    pub m: f64,
249    pub h: f64,
250    pub n: f64,
251    pub a: f64,
252    pub b: f64,
253    pub g_na: f64,
254    pub g_k: f64,
255    pub g_a: f64,
256    pub g_l: f64,
257    pub e_na: f64,
258    pub e_k: f64,
259    pub e_a: f64,
260    pub e_l: f64,
261    pub c_m: f64,
262    pub dt: f64,
263    pub v_threshold: f64,
264}
265
266impl ConnorStevensNeuron {
267    pub fn new() -> Self {
268        Self {
269            v: -68.0,
270            m: 0.01,
271            h: 0.99,
272            n: 0.1,
273            a: 0.5,
274            b: 0.1,
275            g_na: 120.0,
276            g_k: 20.0,
277            g_a: 47.7,
278            g_l: 0.3,
279            e_na: 55.0,
280            e_k: -72.0,
281            e_a: -75.0,
282            e_l: -17.0,
283            c_m: 1.0,
284            dt: 0.01,
285            v_threshold: 0.0,
286        }
287    }
288    pub fn step(&mut self, current: f64) -> i32 {
289        let v_prev = self.v;
290        for _ in 0..10 {
291            let am = safe_rate(0.1, 29.7, self.v, 10.0, 1.0);
292            let bm = 4.0 * (-(self.v + 54.7) / 18.0).exp();
293            let ah = 0.07 * (-(self.v + 48.0) / 20.0).exp();
294            let bh = 1.0 / (1.0 + (-(self.v + 18.0) / 10.0).exp());
295            let an = safe_rate(0.01, 45.7, self.v, 10.0, 0.1);
296            let bn = 0.125 * (-(self.v + 55.7) / 80.0).exp();
297            let a_inf = (0.0761 * ((self.v + 94.22) / 31.84).exp()
298                / (1.0 + ((self.v + 1.17) / 28.93).exp()))
299            .powf(1.0 / 3.0);
300            let b_inf = 1.0 / (1.0 + ((self.v + 53.3) / 14.54).exp()).powf(4.0);
301            let tau_a = 0.3632 + 1.158 / (1.0 + ((self.v + 55.96) / 20.12).exp());
302            let tau_b = 1.24 + 2.678 / (1.0 + ((self.v + 50.0) / 16.027).exp());
303            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
304            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
305            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
306            self.a += (a_inf - self.a) / tau_a * self.dt;
307            self.b += (b_inf - self.b) / tau_b * self.dt;
308            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
309            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
310            let i_a = self.g_a * self.a.powi(3) * self.b * (self.v - self.e_a);
311            let i_l = self.g_l * (self.v - self.e_l);
312            self.v += (-i_na - i_k - i_a - i_l + current) / self.c_m * self.dt;
313        }
314        if self.v >= self.v_threshold && v_prev < self.v_threshold {
315            1
316        } else {
317            0
318        }
319    }
320    pub fn reset(&mut self) {
321        self.v = -68.0;
322        self.m = 0.01;
323        self.h = 0.99;
324        self.n = 0.1;
325        self.a = 0.5;
326        self.b = 0.1;
327    }
328}
329impl Default for ConnorStevensNeuron {
330    fn default() -> Self {
331        Self::new()
332    }
333}
334
335/// Destexhe thalamocortical relay neuron with T-current. Destexhe et al. 1993.
336#[derive(Clone, Debug)]
337pub struct DestexheThalamicNeuron {
338    pub v: f64,
339    pub h_na: f64,
340    pub n_k: f64,
341    pub m_t: f64,
342    pub h_t: f64,
343    pub g_na: f64,
344    pub g_k: f64,
345    pub g_t: f64,
346    pub g_l: f64,
347    pub e_na: f64,
348    pub e_k: f64,
349    pub e_ca: f64,
350    pub e_l: f64,
351    pub dt: f64,
352    pub v_threshold: f64,
353}
354
355impl DestexheThalamicNeuron {
356    pub fn new() -> Self {
357        Self {
358            v: -65.0,
359            h_na: 0.6,
360            n_k: 0.3,
361            m_t: 0.0,
362            h_t: 1.0,
363            g_na: 100.0,
364            g_k: 10.0,
365            g_t: 2.0,
366            g_l: 0.05,
367            e_na: 50.0,
368            e_k: -100.0,
369            e_ca: 120.0,
370            e_l: -70.0,
371            dt: 0.02,
372            v_threshold: -20.0,
373        }
374    }
375    pub fn step(&mut self, current: f64) -> i32 {
376        let v_prev = self.v;
377        for _ in 0..5 {
378            let m_na = 1.0 / (1.0 + (-(self.v + 37.0) / 7.0).exp());
379            let h_na_inf = 1.0 / (1.0 + ((self.v + 41.0) / 4.0).exp());
380            let n_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 12.0).exp());
381            let m_t_inf = 1.0 / (1.0 + (-(self.v + 57.0) / 6.2).exp());
382            let h_t_inf = 1.0 / (1.0 + ((self.v + 81.0) / 4.0).exp());
383            self.h_na += (h_na_inf - self.h_na) / 1.0 * self.dt;
384            self.n_k += (n_inf - self.n_k) / 5.0 * self.dt;
385            self.m_t += (m_t_inf - self.m_t) / 1.0 * self.dt;
386            self.h_t += (h_t_inf - self.h_t) / 20.0 * self.dt;
387            let i_na = self.g_na * m_na.powi(3) * self.h_na * (self.v - self.e_na);
388            let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
389            let i_t = self.g_t * self.m_t.powi(2) * self.h_t * (self.v - self.e_ca);
390            let i_l = self.g_l * (self.v - self.e_l);
391            self.v += (-i_na - i_k - i_t - i_l + current) * self.dt;
392        }
393        if self.v >= self.v_threshold && v_prev < self.v_threshold {
394            1
395        } else {
396            0
397        }
398    }
399    pub fn reset(&mut self) {
400        self.v = -65.0;
401        self.h_na = 0.6;
402        self.n_k = 0.3;
403        self.m_t = 0.0;
404        self.h_t = 1.0;
405    }
406}
407impl Default for DestexheThalamicNeuron {
408    fn default() -> Self {
409        Self::new()
410    }
411}
412
413/// Huber-Braun — temperature-sensitive cold receptor. Braun et al. 1998.
414#[derive(Clone, Debug)]
415pub struct HuberBraunNeuron {
416    pub v: f64,
417    pub a_sd: f64,
418    pub a_sr: f64,
419    pub g_sd: f64,
420    pub g_sr: f64,
421    pub g_l: f64,
422    pub e_sd: f64,
423    pub e_sr: f64,
424    pub e_l: f64,
425    pub tau_sd: f64,
426    pub tau_sr: f64,
427    pub eta: f64,
428    pub dt: f64,
429    pub v_threshold: f64,
430}
431
432impl HuberBraunNeuron {
433    pub fn new() -> Self {
434        Self {
435            v: -50.0,
436            a_sd: 0.0,
437            a_sr: 0.0,
438            g_sd: 1.5,
439            g_sr: 0.4,
440            g_l: 0.1,
441            e_sd: 50.0,
442            e_sr: -90.0,
443            e_l: -60.0,
444            tau_sd: 10.0,
445            tau_sr: 20.0,
446            eta: 0.012,
447            dt: 0.1,
448            v_threshold: -20.0,
449        }
450    }
451    pub fn step(&mut self, current: f64) -> i32 {
452        let v_prev = self.v;
453        let sd_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 5.0).exp());
454        let sr_inf = 1.0 / (1.0 + ((self.v + 40.0) / 5.0).exp());
455        self.a_sd += (sd_inf - self.a_sd) / self.tau_sd * self.dt;
456        self.a_sr += (sr_inf - self.a_sr) / self.tau_sr * self.dt;
457        let i_sd = self.g_sd * self.a_sd * (self.v - self.e_sd);
458        let i_sr = self.g_sr * self.a_sr * (self.v - self.e_sr);
459        let i_l = self.g_l * (self.v - self.e_l);
460        self.v += (-i_sd - i_sr - i_l + current) * self.dt;
461        if self.v >= self.v_threshold && v_prev < self.v_threshold {
462            1
463        } else {
464            0
465        }
466    }
467    pub fn reset(&mut self) {
468        self.v = -50.0;
469        self.a_sd = 0.0;
470        self.a_sr = 0.0;
471    }
472}
473impl Default for HuberBraunNeuron {
474    fn default() -> Self {
475        Self::new()
476    }
477}
478
479/// Golomb fast-spiking interneuron with Kv3. Golomb et al. 2007.
480#[derive(Clone, Debug)]
481pub struct GolombFSNeuron {
482    pub v: f64,
483    pub h: f64,
484    pub n: f64,
485    pub p: f64,
486    pub g_na: f64,
487    pub g_k: f64,
488    pub g_kv3: f64,
489    pub g_l: f64,
490    pub e_na: f64,
491    pub e_k: f64,
492    pub e_l: f64,
493    pub dt: f64,
494    pub v_threshold: f64,
495}
496
497impl GolombFSNeuron {
498    pub fn new() -> Self {
499        Self {
500            v: -65.0,
501            h: 0.9,
502            n: 0.1,
503            p: 0.0,
504            g_na: 112.5,
505            g_k: 225.0,
506            g_kv3: 150.0,
507            g_l: 0.25,
508            e_na: 50.0,
509            e_k: -90.0,
510            e_l: -70.0,
511            dt: 0.01,
512            v_threshold: -20.0,
513        }
514    }
515    pub fn step(&mut self, current: f64) -> i32 {
516        let v_prev = self.v;
517        for _ in 0..10 {
518            let m_inf = 1.0 / (1.0 + (-(self.v + 24.0) / 11.5).exp());
519            let h_inf = 1.0 / (1.0 + ((self.v + 58.3) / 6.7).exp());
520            let n_inf = 1.0 / (1.0 + (-(self.v + 12.4) / 6.8).exp());
521            let p_inf = 1.0 / (1.0 + (-(self.v + 3.0) / 10.0).exp());
522            self.h += (h_inf - self.h) / 0.5 * self.dt;
523            self.n += (n_inf - self.n) / 2.0 * self.dt;
524            self.p += (p_inf - self.p) / 1.0 * self.dt;
525            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.v - self.e_na);
526            let i_k = self.g_k * self.n.powi(2) * (self.v - self.e_k);
527            let i_kv3 = self.g_kv3 * self.p * (self.v - self.e_k);
528            let i_l = self.g_l * (self.v - self.e_l);
529            self.v += (-i_na - i_k - i_kv3 - i_l + current) * self.dt;
530        }
531        if self.v >= self.v_threshold && v_prev < self.v_threshold {
532            1
533        } else {
534            0
535        }
536    }
537    pub fn reset(&mut self) {
538        self.v = -65.0;
539        self.h = 0.9;
540        self.n = 0.1;
541        self.p = 0.0;
542    }
543}
544impl Default for GolombFSNeuron {
545    fn default() -> Self {
546        Self::new()
547    }
548}
549
550/// Pospischil — minimal HH for 5 cortical cell types. Pospischil et al. 2008.
551#[derive(Clone, Debug)]
552pub struct PospischilNeuron {
553    pub v: f64,
554    pub m: f64,
555    pub h: f64,
556    pub n: f64,
557    pub p: f64,
558    pub g_na: f64,
559    pub g_k: f64,
560    pub g_m: f64,
561    pub g_l: f64,
562    pub e_na: f64,
563    pub e_k: f64,
564    pub e_l: f64,
565    pub vt: f64,
566    pub dt: f64,
567    pub v_threshold: f64,
568}
569
570impl PospischilNeuron {
571    pub fn new() -> Self {
572        Self {
573            v: -70.0,
574            m: 0.05,
575            h: 0.6,
576            n: 0.3,
577            p: 0.0,
578            g_na: 50.0,
579            g_k: 5.0,
580            g_m: 0.004,
581            g_l: 0.01,
582            e_na: 50.0,
583            e_k: -90.0,
584            e_l: -70.0,
585            vt: -56.2,
586            dt: 0.025,
587            v_threshold: -20.0,
588        }
589    }
590    pub fn step(&mut self, current: f64) -> i32 {
591        let v_prev = self.v;
592        for _ in 0..4 {
593            let dv = self.v - self.vt;
594            let am = safe_rate(-0.32, 0.0, dv - 13.0, -4.0, 8.0);
595            let bm = safe_rate(0.28, 0.0, dv - 40.0, 5.0, 5.6);
596            let ah = 0.128 * (-(dv - 17.0) / 18.0).exp();
597            let bh = 4.0 / (1.0 + (-(dv - 40.0) / 5.0).exp());
598            let an = safe_rate(-0.032, 0.0, dv - 15.0, -5.0, 0.32);
599            let bn = 0.5 * (-(dv - 10.0) / 40.0).exp();
600            let p_inf = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
601            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
602            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
603            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
604            self.p += (p_inf - self.p) / 100.0 * self.dt;
605            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
606            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
607            let i_m = self.g_m * self.p * (self.v - self.e_k);
608            let i_l = self.g_l * (self.v - self.e_l);
609            self.v += (-i_na - i_k - i_m - i_l + current) * self.dt;
610        }
611        if self.v >= self.v_threshold && v_prev < self.v_threshold {
612            1
613        } else {
614            0
615        }
616    }
617    pub fn reset(&mut self) {
618        self.v = -70.0;
619        self.m = 0.05;
620        self.h = 0.6;
621        self.n = 0.3;
622        self.p = 0.0;
623    }
624}
625impl Default for PospischilNeuron {
626    fn default() -> Self {
627        Self::new()
628    }
629}
630
631/// Mainen-Sejnowski — two-compartment (soma + axon). Mainen & Sejnowski 1996.
632#[derive(Clone, Debug)]
633pub struct MainenSejnowskiNeuron {
634    pub vs: f64,
635    pub va: f64,
636    pub m: f64,
637    pub h: f64,
638    pub n: f64,
639    pub kappa: f64,
640    pub g_na: f64,
641    pub g_k: f64,
642    pub g_l: f64,
643    pub e_na: f64,
644    pub e_k: f64,
645    pub e_l: f64,
646    pub c_s: f64,
647    pub c_a: f64,
648    pub dt: f64,
649    pub v_threshold: f64,
650}
651
652impl MainenSejnowskiNeuron {
653    pub fn new() -> Self {
654        Self {
655            vs: -65.0,
656            va: -65.0,
657            m: 0.05,
658            h: 0.6,
659            n: 0.3,
660            kappa: 10.0,
661            g_na: 3000.0,
662            g_k: 1500.0,
663            g_l: 1.0,
664            e_na: 50.0,
665            e_k: -90.0,
666            e_l: -70.0,
667            c_s: 1.0,
668            c_a: 0.1,
669            dt: 0.005,
670            v_threshold: -20.0,
671        }
672    }
673    pub fn step(&mut self, current: f64) -> i32 {
674        let v_prev = self.vs;
675        for _ in 0..20 {
676            let am = safe_rate(0.1, 40.0, self.va, 10.0, 1.0);
677            let bm = 4.0 * (-(self.va + 65.0) / 18.0).exp();
678            let ah = 0.07 * (-(self.va + 65.0) / 20.0).exp();
679            let bh = 1.0 / (1.0 + (-(self.va + 35.0) / 10.0).exp());
680            let an = safe_rate(0.01, 55.0, self.va, 10.0, 0.1);
681            let bn = 0.125 * (-(self.va + 65.0) / 80.0).exp();
682            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
683            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
684            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
685            let i_na = self.g_na * self.m.powi(3) * self.h * (self.va - self.e_na);
686            let i_k = self.g_k * self.n.powi(4) * (self.va - self.e_k);
687            let i_l_s = self.g_l * (self.vs - self.e_l);
688            let i_c = self.kappa * (self.va - self.vs);
689            self.vs += (-i_l_s + i_c + current) / self.c_s * self.dt;
690            self.va += (-i_na - i_k - self.kappa * (self.va - self.vs)) / self.c_a * self.dt;
691        }
692        if self.vs >= self.v_threshold && v_prev < self.v_threshold {
693            1
694        } else {
695            0
696        }
697    }
698    pub fn reset(&mut self) {
699        self.vs = -65.0;
700        self.va = -65.0;
701        self.m = 0.05;
702        self.h = 0.6;
703        self.n = 0.3;
704    }
705}
706impl Default for MainenSejnowskiNeuron {
707    fn default() -> Self {
708        Self::new()
709    }
710}
711
712/// De Schutter-Bower Purkinje cell — Ca-dependent K. De Schutter & Bower 1994.
713#[derive(Clone, Debug)]
714pub struct DeSchutterPurkinjeNeuron {
715    pub v: f64,
716    pub h_na: f64,
717    pub n_k: f64,
718    pub m_cap: f64,
719    pub h_cap: f64,
720    pub q_kca: f64,
721    pub ca: f64,
722    pub g_na: f64,
723    pub g_k: f64,
724    pub g_cap: f64,
725    pub g_kca: f64,
726    pub g_l: f64,
727    pub e_na: f64,
728    pub e_k: f64,
729    pub e_ca: f64,
730    pub e_l: f64,
731    pub ca_decay: f64,
732    pub f_ca: f64,
733    pub dt: f64,
734    pub v_threshold: f64,
735}
736
737impl DeSchutterPurkinjeNeuron {
738    pub fn new() -> Self {
739        Self {
740            v: -68.0,
741            h_na: 0.8,
742            n_k: 0.1,
743            m_cap: 0.0,
744            h_cap: 0.9,
745            q_kca: 0.0,
746            ca: 0.0001,
747            g_na: 125.0,
748            g_k: 10.0,
749            g_cap: 45.0,
750            g_kca: 35.0,
751            g_l: 0.5,
752            e_na: 45.0,
753            e_k: -85.0,
754            e_ca: 135.0,
755            e_l: -68.0,
756            ca_decay: 0.02,
757            f_ca: 0.00024,
758            dt: 0.01,
759            v_threshold: -20.0,
760        }
761    }
762    pub fn step(&mut self, current: f64) -> i32 {
763        let v_prev = self.v;
764        for _ in 0..5 {
765            let m_na = 1.0 / (1.0 + (-(self.v + 30.0) / 9.0).exp());
766            let h_na_inf = 1.0 / (1.0 + ((self.v + 45.0) / 7.0).exp());
767            let n_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 12.0).exp());
768            let m_cap_inf = 1.0 / (1.0 + (-(self.v + 19.0) / 10.0).exp());
769            let h_cap_inf = 1.0 / (1.0 + ((self.v + 39.0) / 7.0).exp());
770            let q_inf = self.ca / (self.ca + 0.001);
771            self.h_na += (h_na_inf - self.h_na) / 1.0 * self.dt;
772            self.n_k += (n_inf - self.n_k) / 3.0 * self.dt;
773            self.m_cap += (m_cap_inf - self.m_cap) / 1.0 * self.dt;
774            self.h_cap += (h_cap_inf - self.h_cap) / 15.0 * self.dt;
775            self.q_kca += (q_inf - self.q_kca) / 5.0 * self.dt;
776            let i_na = self.g_na * m_na.powi(3) * self.h_na * (self.v - self.e_na);
777            let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
778            let i_cap = self.g_cap * self.m_cap.powi(2) * self.h_cap * (self.v - self.e_ca);
779            let i_kca = self.g_kca * self.q_kca * (self.v - self.e_k);
780            let i_l = self.g_l * (self.v - self.e_l);
781            self.v += (-i_na - i_k - i_cap - i_kca - i_l + current) * self.dt;
782            self.ca = (self.ca + (-self.f_ca * i_cap - self.ca_decay * self.ca) * self.dt).max(0.0);
783        }
784        if self.v >= self.v_threshold && v_prev < self.v_threshold {
785            1
786        } else {
787            0
788        }
789    }
790    pub fn reset(&mut self) {
791        self.v = -68.0;
792        self.h_na = 0.8;
793        self.n_k = 0.1;
794        self.m_cap = 0.0;
795        self.h_cap = 0.9;
796        self.q_kca = 0.0;
797        self.ca = 0.0001;
798    }
799}
800impl Default for DeSchutterPurkinjeNeuron {
801    fn default() -> Self {
802        Self::new()
803    }
804}
805
806/// Plant R15 — Aplysia parabolic burster. Plant & Kim 1976.
807#[derive(Clone, Debug)]
808pub struct PlantR15Neuron {
809    pub v: f64,
810    pub m: f64,
811    pub h: f64,
812    pub n: f64,
813    pub ca: f64,
814    pub g_na: f64,
815    pub g_k: f64,
816    pub g_ca: f64,
817    pub g_l: f64,
818    pub g_kca: f64,
819    pub e_na: f64,
820    pub e_k: f64,
821    pub e_ca: f64,
822    pub e_l: f64,
823    pub c_m: f64,
824    pub k_ca: f64,
825    pub tau_ca: f64,
826    pub dt: f64,
827    pub v_threshold: f64,
828}
829
830impl PlantR15Neuron {
831    pub fn new() -> Self {
832        Self {
833            v: -50.0,
834            m: 0.05,
835            h: 0.6,
836            n: 0.3,
837            ca: 0.1,
838            g_na: 4.0,
839            g_k: 0.3,
840            g_ca: 0.004,
841            g_l: 0.003,
842            g_kca: 0.03,
843            e_na: 30.0,
844            e_k: -75.0,
845            e_ca: 140.0,
846            e_l: -40.0,
847            c_m: 1.0,
848            k_ca: 0.0085,
849            tau_ca: 500.0,
850            dt: 0.05,
851            v_threshold: -10.0,
852        }
853    }
854    pub fn step(&mut self, current: f64) -> i32 {
855        let v_prev = self.v;
856        for _ in 0..5 {
857            let am = safe_rate(0.1, 40.0, self.v, 10.0, 1.0);
858            let bm = 4.0 * (-(self.v + 65.0) / 18.0).exp();
859            let ah = 0.07 * (-(self.v + 65.0) / 20.0).exp();
860            let bh = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
861            let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
862            let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
863            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
864            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
865            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
866            let m_ca = 1.0 / (1.0 + (-(self.v + 25.0) / 5.0).exp());
867            let kca_act = (self.ca / (self.ca + 0.5)).min(1.0);
868            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
869            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
870            let i_ca = self.g_ca * m_ca * (self.v - self.e_ca);
871            let i_kca = self.g_kca * kca_act * (self.v - self.e_k);
872            let i_l = self.g_l * (self.v - self.e_l);
873            self.v += (-i_na - i_k - i_ca - i_kca - i_l + current) / self.c_m * self.dt;
874            self.ca = (self.ca + (-self.k_ca * i_ca - self.ca / self.tau_ca) * self.dt).max(0.0);
875        }
876        if self.v >= self.v_threshold && v_prev < self.v_threshold {
877            1
878        } else {
879            0
880        }
881    }
882    pub fn reset(&mut self) {
883        self.v = -50.0;
884        self.m = 0.05;
885        self.h = 0.6;
886        self.n = 0.3;
887        self.ca = 0.1;
888    }
889}
890impl Default for PlantR15Neuron {
891    fn default() -> Self {
892        Self::new()
893    }
894}
895
896/// Prescott 2008 — Type I/II/III excitability tuning via M-current.
897#[derive(Clone, Debug)]
898pub struct PrescottNeuron {
899    pub v: f64,
900    pub w: f64,
901    pub g_fast: f64,
902    pub g_slow: f64,
903    pub g_l: f64,
904    pub e_fast: f64,
905    pub e_slow: f64,
906    pub e_l: f64,
907    pub beta_w: f64,
908    pub gamma_w: f64,
909    pub tau_w: f64,
910    pub phi: f64,
911    pub dt: f64,
912    pub v_threshold: f64,
913}
914
915impl PrescottNeuron {
916    pub fn new() -> Self {
917        Self {
918            v: -65.0,
919            w: 0.0,
920            g_fast: 20.0,
921            g_slow: 20.0,
922            g_l: 2.0,
923            e_fast: 50.0,
924            e_slow: -100.0,
925            e_l: -70.0,
926            beta_w: -21.0,
927            gamma_w: 15.0,
928            tau_w: 100.0,
929            phi: 0.15,
930            dt: 0.1,
931            v_threshold: -20.0,
932        }
933    }
934    pub fn step(&mut self, current: f64) -> i32 {
935        let v_prev = self.v;
936        let m_inf = 1.0 / (1.0 + (-(self.v + 20.0) / 15.0).exp());
937        let w_inf = 1.0 / (1.0 + (-(self.v - self.beta_w) / self.gamma_w).exp());
938        let i_fast = self.g_fast * m_inf * (self.v - self.e_fast);
939        let i_slow = self.g_slow * self.w * (self.v - self.e_slow);
940        let i_l = self.g_l * (self.v - self.e_l);
941        self.v += (-i_fast - i_slow - i_l + current) * self.dt;
942        self.w += self.phi * (w_inf - self.w) / self.tau_w * self.dt;
943        if self.v >= self.v_threshold && v_prev < self.v_threshold {
944            1
945        } else {
946            0
947        }
948    }
949    pub fn reset(&mut self) {
950        self.v = -65.0;
951        self.w = 0.0;
952    }
953}
954impl Default for PrescottNeuron {
955    fn default() -> Self {
956        Self::new()
957    }
958}
959
960/// Mihalas-Niebur 2009 — generalized IF capturing 20 spike patterns.
961#[derive(Clone, Debug)]
962pub struct MihalasNieburNeuron {
963    pub v: f64,
964    pub theta: f64,
965    pub i1: f64,
966    pub i2: f64,
967    pub v_rest: f64,
968    pub v_reset: f64,
969    pub theta_reset: f64,
970    pub theta_inf: f64,
971    pub tau_v: f64,
972    pub tau_theta: f64,
973    pub tau_1: f64,
974    pub tau_2: f64,
975    pub a: f64,
976    pub b: f64,
977    pub r1: f64,
978    pub r2: f64,
979    pub dt: f64,
980}
981
982impl MihalasNieburNeuron {
983    pub fn new() -> Self {
984        Self {
985            v: 0.0,
986            theta: 1.0,
987            i1: 0.0,
988            i2: 0.0,
989            v_rest: 0.0,
990            v_reset: 0.0,
991            theta_reset: 1.0,
992            theta_inf: 1.0,
993            tau_v: 10.0,
994            tau_theta: 100.0,
995            tau_1: 10.0,
996            tau_2: 200.0,
997            a: 0.0,
998            b: 0.0,
999            r1: 0.0,
1000            r2: 0.0,
1001            dt: 1.0,
1002        }
1003    }
1004    pub fn step(&mut self, current: f64) -> i32 {
1005        self.v += (-(self.v - self.v_rest) + self.i1 + self.i2 + current) / self.tau_v * self.dt;
1006        self.theta += self.a * (self.v - self.v_rest)
1007            + (self.theta_inf - self.theta) / self.tau_theta * self.dt;
1008        self.i1 += -self.i1 / self.tau_1 * self.dt;
1009        self.i2 += -self.i2 / self.tau_2 * self.dt;
1010        if self.v >= self.theta {
1011            self.v = self.v_reset + self.b * (self.v - self.v_rest);
1012            self.theta = self.theta_reset.max(self.theta);
1013            self.i1 += self.r1;
1014            self.i2 += self.r2;
1015            1
1016        } else {
1017            0
1018        }
1019    }
1020    pub fn reset(&mut self) {
1021        self.v = self.v_rest;
1022        self.theta = self.theta_reset;
1023        self.i1 = 0.0;
1024        self.i2 = 0.0;
1025    }
1026}
1027impl Default for MihalasNieburNeuron {
1028    fn default() -> Self {
1029        Self::new()
1030    }
1031}
1032
1033/// Allen GLIF5 — LIF + threshold adaptation + after-spike currents.
1034#[derive(Clone, Debug)]
1035pub struct GLIFNeuron {
1036    pub v: f64,
1037    pub theta: f64,
1038    pub theta_inf: f64,
1039    pub i_asc1: f64,
1040    pub i_asc2: f64,
1041    pub v_rest: f64,
1042    pub v_reset: f64,
1043    pub tau_m: f64,
1044    pub tau_theta: f64,
1045    pub tau_asc1: f64,
1046    pub tau_asc2: f64,
1047    pub a_theta: f64,
1048    pub delta_theta: f64,
1049    pub r_asc1: f64,
1050    pub r_asc2: f64,
1051    pub resistance: f64,
1052    pub dt: f64,
1053}
1054
1055impl GLIFNeuron {
1056    pub fn new() -> Self {
1057        Self {
1058            v: -70.0,
1059            theta: -50.0,
1060            theta_inf: -50.0,
1061            i_asc1: 0.0,
1062            i_asc2: 0.0,
1063            v_rest: -70.0,
1064            v_reset: -70.0,
1065            tau_m: 10.0,
1066            tau_theta: 100.0,
1067            tau_asc1: 10.0,
1068            tau_asc2: 200.0,
1069            a_theta: 0.01,
1070            delta_theta: 2.0,
1071            r_asc1: 1.0,
1072            r_asc2: 0.5,
1073            resistance: 1.0,
1074            dt: 1.0,
1075        }
1076    }
1077    pub fn step(&mut self, current: f64) -> i32 {
1078        self.v += (-(self.v - self.v_rest) + self.resistance * current + self.i_asc1 + self.i_asc2)
1079            / self.tau_m
1080            * self.dt;
1081        self.theta += self.a_theta * (self.v - self.v_rest)
1082            + (self.theta_inf - self.theta) / self.tau_theta * self.dt;
1083        self.i_asc1 *= (-self.dt / self.tau_asc1).exp();
1084        self.i_asc2 *= (-self.dt / self.tau_asc2).exp();
1085        if self.v >= self.theta {
1086            self.v = self.v_reset;
1087            self.theta += self.delta_theta;
1088            self.i_asc1 += self.r_asc1;
1089            self.i_asc2 += self.r_asc2;
1090            1
1091        } else {
1092            0
1093        }
1094    }
1095    pub fn reset(&mut self) {
1096        self.v = self.v_rest;
1097        self.theta = self.theta_inf;
1098        self.i_asc1 = 0.0;
1099        self.i_asc2 = 0.0;
1100    }
1101}
1102impl Default for GLIFNeuron {
1103    fn default() -> Self {
1104        Self::new()
1105    }
1106}
1107
1108/// GIF population — escape-rate generalized IF. Mensi et al. 2012.
1109#[derive(Clone, Debug)]
1110pub struct GIFPopulationNeuron {
1111    pub v: f64,
1112    pub theta: f64,
1113    pub eta: f64,
1114    pub tau_m: f64,
1115    pub tau_eta: f64,
1116    pub delta_v: f64,
1117    pub lambda_0: f64,
1118    pub eta_increment: f64,
1119    pub v_rest: f64,
1120    pub v_reset: f64,
1121    pub dt: f64,
1122    rng: Xoshiro256PlusPlus,
1123}
1124
1125impl GIFPopulationNeuron {
1126    pub fn new(seed: u64) -> Self {
1127        Self {
1128            v: -65.0,
1129            theta: -50.0,
1130            eta: 0.0,
1131            tau_m: 20.0,
1132            tau_eta: 100.0,
1133            delta_v: 2.0,
1134            lambda_0: 0.001,
1135            eta_increment: 5.0,
1136            v_rest: -65.0,
1137            v_reset: -65.0,
1138            dt: 0.5,
1139            rng: Xoshiro256PlusPlus::seed_from_u64(seed),
1140        }
1141    }
1142    pub fn step(&mut self, current: f64) -> i32 {
1143        self.v += (-(self.v - self.v_rest) + current) / self.tau_m * self.dt;
1144        self.eta *= (-self.dt / self.tau_eta).exp();
1145        let hazard = self.lambda_0 * ((self.v - self.theta + self.eta) / self.delta_v).exp();
1146        if self.rng.random::<f64>() < hazard * self.dt {
1147            self.v = self.v_reset;
1148            self.eta += self.eta_increment;
1149            1
1150        } else {
1151            0
1152        }
1153    }
1154    pub fn reset(&mut self) {
1155        self.v = -65.0;
1156        self.eta = 0.0;
1157    }
1158}
1159
1160/// Av-Ron cardiac ganglion — Type III burster. Av-Ron et al. 1991.
1161#[derive(Clone, Debug)]
1162pub struct AvRonCardiacNeuron {
1163    pub v: f64,
1164    pub h: f64,
1165    pub n: f64,
1166    pub s: f64,
1167    pub g_na: f64,
1168    pub g_k: f64,
1169    pub g_s: f64,
1170    pub g_l: f64,
1171    pub e_na: f64,
1172    pub e_k: f64,
1173    pub e_s: f64,
1174    pub e_l: f64,
1175    pub dt: f64,
1176    pub v_threshold: f64,
1177}
1178
1179impl AvRonCardiacNeuron {
1180    pub fn new() -> Self {
1181        Self {
1182            v: -60.0,
1183            h: 0.6,
1184            n: 0.3,
1185            s: 0.5,
1186            g_na: 80.0,
1187            g_k: 40.0,
1188            g_s: 20.0,
1189            g_l: 0.1,
1190            e_na: 40.0,
1191            e_k: -80.0,
1192            e_s: -25.0,
1193            e_l: -60.0,
1194            dt: 0.02,
1195            v_threshold: -20.0,
1196        }
1197    }
1198    pub fn step(&mut self, current: f64) -> i32 {
1199        let v_prev = self.v;
1200        let m_inf = 1.0 / (1.0 + (-(self.v + 35.0) / 7.8).exp());
1201        let h_inf = 1.0 / (1.0 + ((self.v + 55.0) / 7.0).exp());
1202        let n_inf = 1.0 / (1.0 + (-(self.v + 28.0) / 15.0).exp());
1203        let s_inf = 1.0 / (1.0 + (-(self.v + 27.0) / 5.0).exp());
1204        self.h += (h_inf - self.h) / 1.5 * self.dt;
1205        self.n += (n_inf - self.n) / 4.0 * self.dt;
1206        self.s += (s_inf - self.s) / 50.0 * self.dt;
1207        let i_na = self.g_na * m_inf.powi(3) * self.h * (self.v - self.e_na);
1208        let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
1209        let i_s = self.g_s * self.s * (self.v - self.e_s);
1210        let i_l = self.g_l * (self.v - self.e_l);
1211        self.v += (-i_na - i_k - i_s - i_l + current) * self.dt;
1212        if self.v >= self.v_threshold && v_prev < self.v_threshold {
1213            1
1214        } else {
1215            0
1216        }
1217    }
1218    pub fn reset(&mut self) {
1219        self.v = -60.0;
1220        self.h = 0.6;
1221        self.n = 0.3;
1222        self.s = 0.5;
1223    }
1224}
1225impl Default for AvRonCardiacNeuron {
1226    fn default() -> Self {
1227        Self::new()
1228    }
1229}
1230
1231/// Durstewitz PFC neuron with D1 dopamine modulation. Durstewitz et al. 2000.
1232#[derive(Clone, Debug)]
1233pub struct DurstewitzDopamineNeuron {
1234    pub v: f64,
1235    pub h_na: f64,
1236    pub n_k: f64,
1237    pub g_na: f64,
1238    pub g_k: f64,
1239    pub g_nmda: f64,
1240    pub g_l: f64,
1241    pub e_na: f64,
1242    pub e_k: f64,
1243    pub e_nmda: f64,
1244    pub e_l: f64,
1245    pub mg: f64,
1246    pub d1_level: f64,
1247    pub g_nmda_scale: f64,
1248    pub g_k_scale: f64,
1249    pub v_shift_na: f64,
1250    pub dt: f64,
1251    pub v_threshold: f64,
1252}
1253
1254impl DurstewitzDopamineNeuron {
1255    pub fn new() -> Self {
1256        Self {
1257            v: -65.0,
1258            h_na: 0.7,
1259            n_k: 0.2,
1260            g_na: 45.0,
1261            g_k: 18.0,
1262            g_nmda: 0.5,
1263            g_l: 0.02,
1264            e_na: 55.0,
1265            e_k: -80.0,
1266            e_nmda: 0.0,
1267            e_l: -65.0,
1268            mg: 1.0,
1269            d1_level: 0.0,
1270            g_nmda_scale: 2.5,
1271            g_k_scale: 1.5,
1272            v_shift_na: -5.0,
1273            dt: 0.05,
1274            v_threshold: -20.0,
1275        }
1276    }
1277    pub fn step(&mut self, current: f64) -> i32 {
1278        let v_prev = self.v;
1279        let shift = self.d1_level * self.v_shift_na;
1280        let m_inf = 1.0 / (1.0 + (-(self.v + 30.0 + shift) / 9.5).exp());
1281        let h_inf = 1.0 / (1.0 + ((self.v + 53.0) / 7.0).exp());
1282        let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 10.0).exp());
1283        self.h_na += (h_inf - self.h_na) / 1.0 * self.dt;
1284        self.n_k += (n_inf - self.n_k) / 4.0 * self.dt;
1285        let g_eff_nmda = self.g_nmda * (1.0 + self.d1_level * (self.g_nmda_scale - 1.0));
1286        let g_eff_k = self.g_k * (1.0 + self.d1_level * (self.g_k_scale - 1.0));
1287        let mg_block = 1.0 / (1.0 + self.mg * (-0.062 * self.v).exp() / 3.57);
1288        let i_na = self.g_na * m_inf.powi(3) * self.h_na * (self.v - self.e_na);
1289        let i_k = g_eff_k * self.n_k.powi(4) * (self.v - self.e_k);
1290        let i_nmda = g_eff_nmda * mg_block * (self.v - self.e_nmda);
1291        let i_l = self.g_l * (self.v - self.e_l);
1292        self.v += (-i_na - i_k - i_nmda - i_l + current) * self.dt;
1293        if self.v >= self.v_threshold && v_prev < self.v_threshold {
1294            1
1295        } else {
1296            0
1297        }
1298    }
1299    pub fn reset(&mut self) {
1300        self.v = -65.0;
1301        self.h_na = 0.7;
1302        self.n_k = 0.2;
1303    }
1304}
1305impl Default for DurstewitzDopamineNeuron {
1306    fn default() -> Self {
1307        Self::new()
1308    }
1309}
1310
1311/// Hill-Tononi 2005 — thalamocortical sleep/wake with Na-dependent K.
1312#[derive(Clone, Debug)]
1313pub struct HillTononiNeuron {
1314    pub v: f64,
1315    pub h_na: f64,
1316    pub n_k: f64,
1317    pub m_h: f64,
1318    pub h_t: f64,
1319    pub na_i: f64,
1320    pub g_na: f64,
1321    pub g_k: f64,
1322    pub g_h: f64,
1323    pub g_t: f64,
1324    pub g_kna: f64,
1325    pub g_l: f64,
1326    pub e_na: f64,
1327    pub e_k: f64,
1328    pub e_h: f64,
1329    pub e_ca: f64,
1330    pub e_l: f64,
1331    pub na_pump_max: f64,
1332    pub na_eq: f64,
1333    pub dt: f64,
1334    pub v_threshold: f64,
1335}
1336
1337impl HillTononiNeuron {
1338    pub fn new() -> Self {
1339        Self {
1340            v: -65.0,
1341            h_na: 0.6,
1342            n_k: 0.3,
1343            m_h: 0.0,
1344            h_t: 0.9,
1345            na_i: 5.0,
1346            g_na: 50.0,
1347            g_k: 5.0,
1348            g_h: 1.0,
1349            g_t: 3.0,
1350            g_kna: 1.33,
1351            g_l: 0.02,
1352            e_na: 50.0,
1353            e_k: -90.0,
1354            e_h: -43.0,
1355            e_ca: 120.0,
1356            e_l: -70.0,
1357            na_pump_max: 20.0,
1358            na_eq: 9.5,
1359            dt: 0.05,
1360            v_threshold: -20.0,
1361        }
1362    }
1363    pub fn step(&mut self, current: f64) -> i32 {
1364        let v_prev = self.v;
1365        let m_na_inf = 1.0 / (1.0 + (-(self.v + 38.0) / 10.0).exp());
1366        let h_na_inf = 1.0 / (1.0 + ((self.v + 43.0) / 6.0).exp());
1367        let n_k_inf = 1.0 / (1.0 + (-(self.v + 27.0) / 11.5).exp());
1368        let m_h_inf = 1.0 / (1.0 + ((self.v + 75.0) / 5.5).exp());
1369        let m_t_inf = 1.0 / (1.0 + (-(self.v + 59.0) / 6.2).exp());
1370        let h_t_inf = 1.0 / (1.0 + ((self.v + 83.0) / 4.0).exp());
1371        let w_kna = 0.37 / (1.0 + (38.7 / self.na_i.max(0.01)).powf(3.5));
1372        let tau_h_na = (1.0 + 10.0 / (1.0 + ((self.v + 40.0) / 10.0).exp())).max(0.1);
1373        let tau_n_k = (5.0 + 47.0 * (-(((self.v + 50.0) / 25.0).powi(2))))
1374            .max(0.1)
1375            .exp();
1376        let tau_m_h = (20.0
1377            + 1000.0 / (((self.v + 71.5) / 14.2).exp() + (-(self.v + 89.0) / 11.6).exp()))
1378        .max(1.0);
1379        let tau_h_t = if self.v < -81.0 {
1380            (30.8 + 211.4 * ((self.v + 115.2) / 5.0).exp() / (1.0 + ((self.v + 86.0) / 3.2).exp()))
1381                .max(0.1)
1382        } else {
1383            10.0
1384        };
1385        self.h_na += (h_na_inf - self.h_na) / tau_h_na * self.dt;
1386        self.n_k += (n_k_inf - self.n_k) / tau_n_k * self.dt;
1387        self.m_h += (m_h_inf - self.m_h) / tau_m_h * self.dt;
1388        self.h_t += (h_t_inf - self.h_t) / tau_h_t * self.dt;
1389        let i_na = self.g_na * m_na_inf.powi(3) * self.h_na * (self.v - self.e_na);
1390        let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
1391        let i_h = self.g_h * self.m_h * (self.v - self.e_h);
1392        let i_t = self.g_t * m_t_inf.powi(2) * self.h_t * (self.v - self.e_ca);
1393        let i_kna = self.g_kna * w_kna * (self.v - self.e_k);
1394        let i_l = self.g_l * (self.v - self.e_l);
1395        self.v += (-i_na - i_k - i_h - i_t - i_kna - i_l + current) * self.dt;
1396        self.na_i = (self.na_i
1397            + (-0.001 * i_na - self.na_pump_max * (self.na_i / (self.na_i + self.na_eq)))
1398                * self.dt)
1399            .max(0.0);
1400        if self.v >= self.v_threshold && v_prev < self.v_threshold {
1401            1
1402        } else {
1403            0
1404        }
1405    }
1406    pub fn reset(&mut self) {
1407        self.v = -65.0;
1408        self.h_na = 0.6;
1409        self.n_k = 0.3;
1410        self.m_h = 0.0;
1411        self.h_t = 0.9;
1412        self.na_i = 5.0;
1413    }
1414}
1415impl Default for HillTononiNeuron {
1416    fn default() -> Self {
1417        Self::new()
1418    }
1419}
1420
1421/// Bertram phantom burster — dual slow K for phantom bursting. Bertram et al. 2000.
1422#[derive(Clone, Debug)]
1423pub struct BertramPhantomBurster {
1424    pub v: f64,
1425    pub s1: f64,
1426    pub s2: f64,
1427    pub g_ca: f64,
1428    pub g_k: f64,
1429    pub g_s1: f64,
1430    pub g_s2: f64,
1431    pub g_l: f64,
1432    pub e_ca: f64,
1433    pub e_k: f64,
1434    pub e_l: f64,
1435    pub c_m: f64,
1436    pub v_m: f64,
1437    pub s_m: f64,
1438    pub v_n: f64,
1439    pub s_n: f64,
1440    pub v_s1: f64,
1441    pub s_s1: f64,
1442    pub v_s2: f64,
1443    pub s_s2: f64,
1444    pub tau_s1: f64,
1445    pub tau_s2: f64,
1446    pub dt: f64,
1447    pub v_threshold: f64,
1448}
1449
1450impl BertramPhantomBurster {
1451    pub fn new() -> Self {
1452        Self {
1453            v: -50.0,
1454            s1: 0.1,
1455            s2: 0.1,
1456            g_ca: 3.6,
1457            g_k: 10.0,
1458            g_s1: 4.0,
1459            g_s2: 4.0,
1460            g_l: 0.2,
1461            e_ca: 25.0,
1462            e_k: -75.0,
1463            e_l: -40.0,
1464            c_m: 5.3,
1465            v_m: -20.0,
1466            s_m: 12.0,
1467            v_n: -16.0,
1468            s_n: 5.6,
1469            v_s1: -35.0,
1470            s_s1: 10.0,
1471            v_s2: -35.0,
1472            s_s2: 10.0,
1473            tau_s1: 20000.0,
1474            tau_s2: 100000.0,
1475            dt: 0.5,
1476            v_threshold: -20.0,
1477        }
1478    }
1479    pub fn step(&mut self, current: f64) -> i32 {
1480        let v_prev = self.v;
1481        let m_inf = 1.0 / (1.0 + (-(self.v - self.v_m) / self.s_m).exp());
1482        let n_inf = 1.0 / (1.0 + (-(self.v - self.v_n) / self.s_n).exp());
1483        let s1_inf = 1.0 / (1.0 + (-(self.v - self.v_s1) / self.s_s1).exp());
1484        let s2_inf = 1.0 / (1.0 + (-(self.v - self.v_s2) / self.s_s2).exp());
1485        let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
1486        let i_k = self.g_k * n_inf * (self.v - self.e_k);
1487        let i_s1 = self.g_s1 * self.s1 * (self.v - self.e_k);
1488        let i_s2 = self.g_s2 * self.s2 * (self.v - self.e_k);
1489        let i_l = self.g_l * (self.v - self.e_l);
1490        self.v += (-i_ca - i_k - i_s1 - i_s2 - i_l + current) / self.c_m * self.dt;
1491        self.s1 += (s1_inf - self.s1) / self.tau_s1 * self.dt;
1492        self.s2 += (s2_inf - self.s2) / self.tau_s2 * self.dt;
1493        if self.v >= self.v_threshold && v_prev < self.v_threshold {
1494            1
1495        } else {
1496            0
1497        }
1498    }
1499    pub fn reset(&mut self) {
1500        self.v = -50.0;
1501        self.s1 = 0.1;
1502        self.s2 = 0.1;
1503    }
1504}
1505impl Default for BertramPhantomBurster {
1506    fn default() -> Self {
1507        Self::new()
1508    }
1509}
1510
1511/// Yamada 1989 — subcritical Hopf burster.
1512#[derive(Clone, Debug)]
1513pub struct YamadaNeuron {
1514    pub v: f64,
1515    pub n: f64,
1516    pub q: f64,
1517    pub g_na: f64,
1518    pub g_k: f64,
1519    pub g_q: f64,
1520    pub g_l: f64,
1521    pub e_na: f64,
1522    pub e_k: f64,
1523    pub e_q: f64,
1524    pub e_l: f64,
1525    pub tau_q: f64,
1526    pub dt: f64,
1527    pub v_threshold: f64,
1528}
1529
1530impl YamadaNeuron {
1531    pub fn new() -> Self {
1532        Self {
1533            v: -60.0,
1534            n: 0.1,
1535            q: 0.0,
1536            g_na: 20.0,
1537            g_k: 10.0,
1538            g_q: 5.0,
1539            g_l: 0.5,
1540            e_na: 60.0,
1541            e_k: -80.0,
1542            e_q: -80.0,
1543            e_l: -60.0,
1544            tau_q: 300.0,
1545            dt: 0.05,
1546            v_threshold: -20.0,
1547        }
1548    }
1549    pub fn step(&mut self, current: f64) -> i32 {
1550        let v_prev = self.v;
1551        let m_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 9.5).exp());
1552        let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 10.0).exp());
1553        let q_inf = 1.0 / (1.0 + (-(self.v + 50.0) / 10.0).exp());
1554        let tau_n = 1.0 + 7.5 / (1.0 + ((self.v + 40.0) / 12.0).exp());
1555        let i_na = self.g_na * m_inf.powi(3) * (1.0 - self.n) * (self.v - self.e_na);
1556        let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
1557        let i_q = self.g_q * self.q * (self.v - self.e_q);
1558        let i_l = self.g_l * (self.v - self.e_l);
1559        self.v += (-i_na - i_k - i_q - i_l + current) * self.dt;
1560        self.n += (n_inf - self.n) / tau_n * self.dt;
1561        self.q += (q_inf - self.q) / self.tau_q * self.dt;
1562        if self.v >= self.v_threshold && v_prev < self.v_threshold {
1563            1
1564        } else {
1565            0
1566        }
1567    }
1568    pub fn reset(&mut self) {
1569        self.v = -60.0;
1570        self.n = 0.1;
1571        self.q = 0.0;
1572    }
1573}
1574impl Default for YamadaNeuron {
1575    fn default() -> Self {
1576        Self::new()
1577    }
1578}
1579
1580#[cfg(test)]
1581mod tests {
1582    use super::*;
1583
1584    #[test]
1585    fn hh_fires() {
1586        let mut n = HodgkinHuxleyNeuron::new();
1587        let t: i32 = (0..100).map(|_| n.step(10.0)).sum();
1588        assert!(t > 0);
1589    }
1590    #[test]
1591    fn traub_fires() {
1592        let mut n = TraubMilesNeuron::new();
1593        let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
1594        assert!(t > 0);
1595    }
1596    #[test]
1597    fn wb_fires() {
1598        let mut n = WangBuzsakiNeuron::new();
1599        let t: i32 = (0..200).map(|_| n.step(2.0)).sum();
1600        assert!(t > 0);
1601    }
1602    #[test]
1603    fn cs_fires() {
1604        let mut n = ConnorStevensNeuron::new();
1605        let t: i32 = (0..200).map(|_| n.step(10.0)).sum();
1606        assert!(t > 0);
1607    }
1608    #[test]
1609    fn destexhe_fires() {
1610        let mut n = DestexheThalamicNeuron::new();
1611        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
1612        assert!(t > 0);
1613    }
1614    #[test]
1615    fn hb_fires() {
1616        let mut n = HuberBraunNeuron::new();
1617        let t: i32 = (0..5000).map(|_| n.step(10.0)).sum();
1618        assert!(t > 0);
1619    }
1620    #[test]
1621    fn golomb_fires() {
1622        let mut n = GolombFSNeuron::new();
1623        let t: i32 = (0..2000).map(|_| n.step(200.0)).sum();
1624        assert!(t > 0);
1625    }
1626    #[test]
1627    fn pospischil_fires() {
1628        let mut n = PospischilNeuron::new();
1629        let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
1630        assert!(t > 0);
1631    }
1632    #[test]
1633    fn mainen_fires() {
1634        let mut n = MainenSejnowskiNeuron::new();
1635        let t: i32 = (0..5000).map(|_| n.step(500.0)).sum();
1636        assert!(t > 0);
1637    }
1638    #[test]
1639    fn purkinje_fires() {
1640        let mut n = DeSchutterPurkinjeNeuron::new();
1641        let t: i32 = (0..5000).map(|_| n.step(50.0)).sum();
1642        assert!(t > 0);
1643    }
1644    #[test]
1645    fn plant_r15_fires() {
1646        let mut n = PlantR15Neuron::new();
1647        let t: i32 = (0..500).map(|_| n.step(2.0)).sum();
1648        assert!(t > 0);
1649    }
1650    #[test]
1651    fn prescott_fires() {
1652        let mut n = PrescottNeuron::new();
1653        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
1654        assert!(t > 0);
1655    }
1656    #[test]
1657    fn mn_fires() {
1658        let mut n = MihalasNieburNeuron::new();
1659        let t: i32 = (0..100).map(|_| n.step(5.0)).sum();
1660        assert!(t > 0);
1661    }
1662    #[test]
1663    fn glif_fires() {
1664        let mut n = GLIFNeuron::new();
1665        let t: i32 = (0..200).map(|_| n.step(30.0)).sum();
1666        assert!(t > 0);
1667    }
1668    #[test]
1669    fn gif_pop_fires() {
1670        let mut n = GIFPopulationNeuron::new(42);
1671        let t: i32 = (0..1000).map(|_| n.step(30.0)).sum();
1672        assert!(t > 0);
1673    }
1674    #[test]
1675    fn avron_fires() {
1676        let mut n = AvRonCardiacNeuron::new();
1677        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
1678        assert!(t > 0);
1679    }
1680    #[test]
1681    fn durstewitz_fires() {
1682        let mut n = DurstewitzDopamineNeuron::new();
1683        let t: i32 = (0..1000).map(|_| n.step(3.0)).sum();
1684        assert!(t > 0);
1685    }
1686    #[test]
1687    fn hill_tononi_fires() {
1688        let mut n = HillTononiNeuron::new();
1689        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
1690        assert!(t > 0);
1691    }
1692    #[test]
1693    fn bertram_fires() {
1694        let mut n = BertramPhantomBurster::new();
1695        let t: i32 = (0..10000).map(|_| n.step(200.0)).sum();
1696        assert!(t > 0);
1697    }
1698    #[test]
1699    fn yamada_fires() {
1700        let mut n = YamadaNeuron::new();
1701        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
1702        assert!(t > 0);
1703    }
1704}
sc_neurocore_engine/neurons/biophysical.rs

sc_neurocore_engine/neurons/
biophysical.rs
