1use rand::{RngExt, SeedableRng};
12use rand_xoshiro::Xoshiro256PlusPlus;
13
14pub fn safe_rate(
17 number_factor: f64,
18 v_offset: f64,
19 v: f64,
20 denom_scale: f64,
21 fallback: f64,
22) -> f64 {
23 let d = v + v_offset;
24 if d.abs() < 1e-7 {
25 fallback
26 } else {
27 number_factor * d / (1.0 - (-d / denom_scale).exp())
28 }
29}
30
31#[derive(Clone, Debug)]
33pub struct HodgkinHuxleyNeuron {
34 pub v: f64,
35 pub m: f64,
36 pub h: f64,
37 pub n: f64,
38 pub c_m: f64,
39 pub g_na: f64,
40 pub g_k: f64,
41 pub g_l: f64,
42 pub e_na: f64,
43 pub e_k: f64,
44 pub e_l: f64,
45 pub dt: f64,
46 pub v_threshold: f64,
47}
48
49impl HodgkinHuxleyNeuron {
50 pub fn new() -> Self {
51 Self {
52 v: -65.0,
53 m: 0.05,
54 h: 0.6,
55 n: 0.32,
56 c_m: 1.0,
57 g_na: 120.0,
58 g_k: 36.0,
59 g_l: 0.3,
60 e_na: 50.0,
61 e_k: -77.0,
62 e_l: -54.4,
63 dt: 0.01,
64 v_threshold: 0.0,
65 }
66 }
67 pub fn step(&mut self, current: f64) -> i32 {
68 let v_prev = self.v;
69 let steps = (1.0 / self.dt) as usize;
70 for _ in 0..steps {
71 let am = safe_rate(0.1, 40.0, self.v, 10.0, 1.0);
72 let bm = 4.0 * (-(self.v + 65.0) / 18.0).exp();
73 let ah = 0.07 * (-(self.v + 65.0) / 20.0).exp();
74 let bh = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
75 let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
76 let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
77 self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
78 self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
79 self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
80 let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
81 let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
82 let i_l = self.g_l * (self.v - self.e_l);
83 self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
84 }
85 if self.v >= self.v_threshold && v_prev < self.v_threshold {
86 1
87 } else {
88 0
89 }
90 }
91 pub fn reset(&mut self) {
92 self.v = -65.0;
93 self.m = 0.05;
94 self.h = 0.6;
95 self.n = 0.32;
96 }
97}
98impl Default for HodgkinHuxleyNeuron {
99 fn default() -> Self {
100 Self::new()
101 }
102}
103
104#[derive(Clone, Debug)]
121pub struct TraubMilesNeuron {
122 pub v: f64,
123 pub m: f64,
124 pub h: f64,
125 pub n: f64,
126 pub g_na: f64,
127 pub g_k: f64,
128 pub g_l: f64,
129 pub e_na: f64,
130 pub e_k: f64,
131 pub e_l: f64,
132 pub dt: f64,
133 pub v_threshold: f64,
134}
135
136impl TraubMilesNeuron {
137 pub fn new() -> Self {
138 Self {
139 v: -67.0,
140 m: 0.05,
141 h: 0.6,
142 n: 0.3,
143 g_na: 100.0,
144 g_k: 80.0,
145 g_l: 0.1,
146 e_na: 50.0,
147 e_k: -100.0,
148 e_l: -67.0,
149 dt: 0.01,
150 v_threshold: -20.0,
151 }
152 }
153 fn finite_gate(value: f64) -> bool {
154 value.is_finite() && (0.0..=1.0).contains(&value)
155 }
156 fn valid_runtime(&self) -> bool {
157 self.v.is_finite()
158 && Self::finite_gate(self.m)
159 && Self::finite_gate(self.h)
160 && Self::finite_gate(self.n)
161 && self.g_na.is_finite()
162 && self.g_na >= 0.0
163 && self.g_k.is_finite()
164 && self.g_k >= 0.0
165 && self.g_l.is_finite()
166 && self.g_l >= 0.0
167 && self.e_na.is_finite()
168 && self.e_k.is_finite()
169 && self.e_l.is_finite()
170 && self.dt.is_finite()
171 && self.dt > 0.0
172 && self.v_threshold.is_finite()
173 }
174 fn rates(v: f64) -> Option<(f64, f64, f64, f64, f64, f64)> {
175 let d = v + 54.0;
176 let am = if d.abs() > 1e-6 {
177 0.32 * d / (1.0 - (-d / 4.0).exp())
178 } else {
179 8.0
180 };
181 let d2 = v + 27.0;
182 let bm = if d2.abs() > 1e-6 {
183 0.28 * d2 / ((d2 / 5.0).exp() - 1.0)
184 } else {
185 5.6
186 };
187 let ah = 0.128 * (-(v + 50.0) / 18.0).exp();
188 let bh = 4.0 / (1.0 + (-(v + 27.0) / 5.0).exp());
189 let d3 = v + 52.0;
190 let an = if d3.abs() > 1e-6 {
191 0.032 * d3 / (1.0 - (-d3 / 5.0).exp())
192 } else {
193 0.32
194 };
195 let bn = 0.5 * (-(v + 57.0) / 40.0).exp();
196 if [am, bm, ah, bh, an, bn]
197 .iter()
198 .all(|rate| rate.is_finite() && *rate >= 0.0)
199 {
200 Some((am, bm, ah, bh, an, bn))
201 } else {
202 None
203 }
204 }
205 fn derivatives(
206 &self,
207 v: f64,
208 m: f64,
209 h: f64,
210 n: f64,
211 current: f64,
212 ) -> Option<(f64, f64, f64, f64)> {
213 if !v.is_finite() || !Self::finite_gate(m) || !Self::finite_gate(h) || !Self::finite_gate(n)
214 {
215 return None;
216 }
217 let (am, bm, ah, bh, an, bn) = Self::rates(v)?;
218 let dm = am * (1.0 - m) - bm * m;
219 let dh = ah * (1.0 - h) - bh * h;
220 let dn = an * (1.0 - n) - bn * n;
221 let i_na = self.g_na * m.powi(3) * h * (v - self.e_na);
222 let i_k = self.g_k * n.powi(4) * (v - self.e_k);
223 let i_l = self.g_l * (v - self.e_l);
224 let dv = -i_na - i_k - i_l + current;
225 if [dv, dm, dh, dn, i_na, i_k, i_l]
226 .iter()
227 .all(|value| value.is_finite())
228 {
229 Some((dv, dm, dh, dn))
230 } else {
231 None
232 }
233 }
234 fn rk4_substep(
235 &self,
236 v: f64,
237 m: f64,
238 h: f64,
239 n: f64,
240 current: f64,
241 ) -> Option<(f64, f64, f64, f64)> {
242 let (k1_v, k1_m, k1_h, k1_n) = self.derivatives(v, m, h, n, current)?;
243 let (k2_v, k2_m, k2_h, k2_n) = self.derivatives(
244 v + 0.5 * self.dt * k1_v,
245 m + 0.5 * self.dt * k1_m,
246 h + 0.5 * self.dt * k1_h,
247 n + 0.5 * self.dt * k1_n,
248 current,
249 )?;
250 let (k3_v, k3_m, k3_h, k3_n) = self.derivatives(
251 v + 0.5 * self.dt * k2_v,
252 m + 0.5 * self.dt * k2_m,
253 h + 0.5 * self.dt * k2_h,
254 n + 0.5 * self.dt * k2_n,
255 current,
256 )?;
257 let (k4_v, k4_m, k4_h, k4_n) = self.derivatives(
258 v + self.dt * k3_v,
259 m + self.dt * k3_m,
260 h + self.dt * k3_h,
261 n + self.dt * k3_n,
262 current,
263 )?;
264 let next_v = v + self.dt * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v) / 6.0;
265 let next_m = m + self.dt * (k1_m + 2.0 * k2_m + 2.0 * k3_m + k4_m) / 6.0;
266 let next_h = h + self.dt * (k1_h + 2.0 * k2_h + 2.0 * k3_h + k4_h) / 6.0;
267 let next_n = n + self.dt * (k1_n + 2.0 * k2_n + 2.0 * k3_n + k4_n) / 6.0;
268 if next_v.is_finite()
269 && Self::finite_gate(next_m)
270 && Self::finite_gate(next_h)
271 && Self::finite_gate(next_n)
272 {
273 Some((next_v, next_m, next_h, next_n))
274 } else {
275 None
276 }
277 }
278 pub fn step(&mut self, current: f64) -> i32 {
279 if !current.is_finite() || !self.valid_runtime() {
280 return 0;
281 }
282 let v_prev = self.v;
283 let mut v = self.v;
284 let mut m = self.m;
285 let mut h = self.h;
286 let mut n = self.n;
287 for _ in 0..10 {
288 let Some((next_v, next_m, next_h, next_n)) = self.rk4_substep(v, m, h, n, current)
289 else {
290 return 0;
291 };
292 v = next_v;
293 m = next_m;
294 h = next_h;
295 n = next_n;
296 }
297 self.v = v;
298 self.m = m;
299 self.h = h;
300 self.n = n;
301 if self.v >= self.v_threshold && v_prev < self.v_threshold {
302 1
303 } else {
304 0
305 }
306 }
307 pub fn reset(&mut self) {
308 *self = Self::new();
309 }
310}
311impl Default for TraubMilesNeuron {
312 fn default() -> Self {
313 Self::new()
314 }
315}
316
317#[derive(Clone, Debug)]
319pub struct WangBuzsakiNeuron {
320 pub v: f64,
321 pub h: f64,
322 pub n: f64,
323 pub g_na: f64,
324 pub g_k: f64,
325 pub g_l: f64,
326 pub e_na: f64,
327 pub e_k: f64,
328 pub e_l: f64,
329 pub c_m: f64,
330 pub phi: f64,
331 pub dt: f64,
332 pub v_threshold: f64,
333}
334
335impl WangBuzsakiNeuron {
336 pub fn new() -> Self {
337 Self {
338 v: -65.0,
339 h: 0.8,
340 n: 0.1,
341 g_na: 35.0,
342 g_k: 9.0,
343 g_l: 0.1,
344 e_na: 55.0,
345 e_k: -90.0,
346 e_l: -65.0,
347 c_m: 1.0,
348 phi: 5.0,
349 dt: 0.01,
350 v_threshold: -20.0,
351 }
352 }
353 pub fn step(&mut self, current: f64) -> i32 {
354 let v_prev = self.v;
355 let n_sub = (0.5 / self.dt.max(0.001)) as usize;
356 for _ in 0..n_sub {
357 let am = safe_rate(0.1, 35.0, self.v, 10.0, 1.0);
358 let bm = 4.0 * (-(self.v + 60.0) / 18.0).exp();
359 let m_inf = am / (am + bm);
360 let ah = 0.07 * (-(self.v + 58.0) / 20.0).exp();
361 let bh = 1.0 / (1.0 + (-(self.v + 28.0) / 10.0).exp());
362 let an = safe_rate(0.01, 34.0, self.v, 10.0, 0.1);
363 let bn = 0.125 * (-(self.v + 44.0) / 80.0).exp();
364 self.h += self.phi * (ah * (1.0 - self.h) - bh * self.h) * self.dt;
365 self.n += self.phi * (an * (1.0 - self.n) - bn * self.n) * self.dt;
366 let i_na = self.g_na * m_inf.powi(3) * self.h * (self.v - self.e_na);
367 let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
368 let i_l = self.g_l * (self.v - self.e_l);
369 self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
370 }
371 if self.v >= self.v_threshold && v_prev < self.v_threshold {
372 1
373 } else {
374 0
375 }
376 }
377 pub fn reset(&mut self) {
378 self.v = -65.0;
379 self.h = 0.8;
380 self.n = 0.1;
381 }
382}
383impl Default for WangBuzsakiNeuron {
384 fn default() -> Self {
385 Self::new()
386 }
387}
388
389#[derive(Clone, Debug)]
391pub struct ConnorStevensNeuron {
392 pub v: f64,
393 pub m: f64,
394 pub h: f64,
395 pub n: f64,
396 pub a: f64,
397 pub b: f64,
398 pub g_na: f64,
399 pub g_k: f64,
400 pub g_a: f64,
401 pub g_l: f64,
402 pub e_na: f64,
403 pub e_k: f64,
404 pub e_a: f64,
405 pub e_l: f64,
406 pub c_m: f64,
407 pub dt: f64,
408 pub v_threshold: f64,
409}
410
411impl ConnorStevensNeuron {
412 pub fn new() -> Self {
413 Self {
414 v: -68.0,
415 m: 0.01,
416 h: 0.99,
417 n: 0.1,
418 a: 0.5,
419 b: 0.1,
420 g_na: 120.0,
421 g_k: 20.0,
422 g_a: 47.7,
423 g_l: 0.3,
424 e_na: 55.0,
425 e_k: -72.0,
426 e_a: -75.0,
427 e_l: -17.0,
428 c_m: 1.0,
429 dt: 0.01,
430 v_threshold: 0.0,
431 }
432 }
433 fn cs_safe_exp(x: f64) -> Option<f64> {
434 if x.is_finite() && x <= 700.0 {
435 Some(x.exp())
436 } else {
437 None
438 }
439 }
440
441 fn cs_safe_rate(scale: f64, shift: f64, v: f64, denom: f64) -> Option<f64> {
442 let delta = v + shift;
443 let x = delta / denom;
444 if x.abs() < 1e-9 {
445 return Some(scale * denom);
446 }
447 let e = Self::cs_safe_exp(-x)?;
448 let value = scale * delta / (1.0 - e);
449 value.is_finite().then_some(value)
450 }
451
452 fn cs_valid_state(v: f64, m: f64, h: f64, n: f64, a: f64, b: f64) -> bool {
453 [v, m, h, n, a, b].iter().all(|x| x.is_finite())
454 && (-250.0..=250.0).contains(&v)
455 && (-0.05..=1.05).contains(&m)
456 && (-0.05..=1.05).contains(&h)
457 && (-0.05..=1.05).contains(&n)
458 && (-0.05..=1.5).contains(&a)
459 && (-0.05..=1.05).contains(&b)
460 }
461
462 fn cs_valid_static(&self) -> bool {
463 [
464 self.g_na,
465 self.g_k,
466 self.g_a,
467 self.g_l,
468 self.e_na,
469 self.e_k,
470 self.e_a,
471 self.e_l,
472 self.c_m,
473 self.dt,
474 self.v_threshold,
475 ]
476 .iter()
477 .all(|x| x.is_finite())
478 && self.g_na >= 0.0
479 && self.g_k >= 0.0
480 && self.g_a >= 0.0
481 && self.g_l >= 0.0
482 && self.c_m > 0.0
483 && self.dt > 0.0
484 }
485
486 fn cs_derivatives(
487 &self,
488 state: (f64, f64, f64, f64, f64, f64),
489 current: f64,
490 ) -> Option<(f64, f64, f64, f64, f64, f64)> {
491 let (v, m, h, n, a, b) = state;
492 let am = Self::cs_safe_rate(0.38, 29.7, v, 10.0)?;
493 let bm = 15.2 * Self::cs_safe_exp(-(v + 54.7) / 18.0)?;
494 let ah = 0.266 * Self::cs_safe_exp(-(v + 48.0) / 20.0)?;
495 let bh = 3.8 / (1.0 + Self::cs_safe_exp(-(v + 18.0) / 10.0)?);
496 let an = Self::cs_safe_rate(0.02, 45.7, v, 10.0)?;
497 let bn = 0.25 * Self::cs_safe_exp(-(v + 55.7) / 80.0)?;
498 let a_base = 0.0761 * Self::cs_safe_exp((v + 94.22) / 31.84)?
499 / (1.0 + Self::cs_safe_exp((v + 1.17) / 28.93)?);
500 if !a_base.is_finite() || a_base < 0.0 {
501 return None;
502 }
503 let a_inf = a_base.powf(1.0 / 3.0);
504 let tau_a = 0.3632 + 1.158 / (1.0 + Self::cs_safe_exp((v + 55.96) / 20.12)?);
505 let b_base = 1.0 / (1.0 + Self::cs_safe_exp((v + 53.3) / 14.54)?);
506 let b_inf = b_base.powf(4.0);
507 let tau_b = 1.24 + 2.678 / (1.0 + Self::cs_safe_exp((v + 50.0) / 16.027)?);
508 let i_na = self.g_na * m.powi(3) * h * (v - self.e_na);
509 let i_k = self.g_k * n.powi(4) * (v - self.e_k);
510 let i_a = self.g_a * a.powi(3) * b * (v - self.e_a);
511 let i_l = self.g_l * (v - self.e_l);
512 let deriv = (
513 (-i_na - i_k - i_a - i_l + current) / self.c_m,
514 am * (1.0 - m) - bm * m,
515 ah * (1.0 - h) - bh * h,
516 an * (1.0 - n) - bn * n,
517 (a_inf - a) / tau_a,
518 (b_inf - b) / tau_b,
519 );
520 [deriv.0, deriv.1, deriv.2, deriv.3, deriv.4, deriv.5]
521 .iter()
522 .all(|x| x.is_finite())
523 .then_some(deriv)
524 }
525
526 fn cs_rk4_candidate(&self, current: f64) -> Option<(f64, f64, f64, f64, f64, f64)> {
527 if !self.cs_valid_static()
528 || !current.is_finite()
529 || !Self::cs_valid_state(self.v, self.m, self.h, self.n, self.a, self.b)
530 {
531 return None;
532 }
533 let mut state = (self.v, self.m, self.h, self.n, self.a, self.b);
534 let substeps = (1.0 / self.dt.max(0.001)) as usize;
535 for _ in 0..substeps {
536 let k1 = self.cs_derivatives(state, current)?;
537 let k2_state = (
538 state.0 + 0.5 * self.dt * k1.0,
539 state.1 + 0.5 * self.dt * k1.1,
540 state.2 + 0.5 * self.dt * k1.2,
541 state.3 + 0.5 * self.dt * k1.3,
542 state.4 + 0.5 * self.dt * k1.4,
543 state.5 + 0.5 * self.dt * k1.5,
544 );
545 let k2 = self.cs_derivatives(k2_state, current)?;
546 let k3_state = (
547 state.0 + 0.5 * self.dt * k2.0,
548 state.1 + 0.5 * self.dt * k2.1,
549 state.2 + 0.5 * self.dt * k2.2,
550 state.3 + 0.5 * self.dt * k2.3,
551 state.4 + 0.5 * self.dt * k2.4,
552 state.5 + 0.5 * self.dt * k2.5,
553 );
554 let k3 = self.cs_derivatives(k3_state, current)?;
555 let k4_state = (
556 state.0 + self.dt * k3.0,
557 state.1 + self.dt * k3.1,
558 state.2 + self.dt * k3.2,
559 state.3 + self.dt * k3.3,
560 state.4 + self.dt * k3.4,
561 state.5 + self.dt * k3.5,
562 );
563 let k4 = self.cs_derivatives(k4_state, current)?;
564 state = (
565 state.0 + self.dt * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) / 6.0,
566 state.1 + self.dt * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1) / 6.0,
567 state.2 + self.dt * (k1.2 + 2.0 * k2.2 + 2.0 * k3.2 + k4.2) / 6.0,
568 state.3 + self.dt * (k1.3 + 2.0 * k2.3 + 2.0 * k3.3 + k4.3) / 6.0,
569 state.4 + self.dt * (k1.4 + 2.0 * k2.4 + 2.0 * k3.4 + k4.4) / 6.0,
570 state.5 + self.dt * (k1.5 + 2.0 * k2.5 + 2.0 * k3.5 + k4.5) / 6.0,
571 );
572 if !Self::cs_valid_state(state.0, state.1, state.2, state.3, state.4, state.5) {
573 return None;
574 }
575 }
576 Some(state)
577 }
578
579 pub fn step(&mut self, current: f64) -> i32 {
580 let v_prev = self.v;
581 let Some((v, m, h, n, a, b)) = self.cs_rk4_candidate(current) else {
582 return 0;
583 };
584 self.v = v;
585 self.m = m;
586 self.h = h;
587 self.n = n;
588 self.a = a;
589 self.b = b;
590 if self.v >= self.v_threshold && v_prev < self.v_threshold {
591 1
592 } else {
593 0
594 }
595 }
596 pub fn reset(&mut self) {
597 self.v = -68.0;
598 self.m = 0.01;
599 self.h = 0.99;
600 self.n = 0.1;
601 self.a = 0.5;
602 self.b = 0.1;
603 }
604}
605impl Default for ConnorStevensNeuron {
606 fn default() -> Self {
607 Self::new()
608 }
609}
610
611#[derive(Clone, Debug)]
613pub struct DestexheThalamicNeuron {
614 pub v: f64,
615 pub h_na: f64,
616 pub n_k: f64,
617 pub m_t: f64,
618 pub h_t: f64,
619 pub g_na: f64,
620 pub g_k: f64,
621 pub g_t: f64,
622 pub g_l: f64,
623 pub e_na: f64,
624 pub e_k: f64,
625 pub e_ca: f64,
626 pub e_l: f64,
627 pub dt: f64,
628 pub v_threshold: f64,
629}
630
631impl DestexheThalamicNeuron {
632 pub fn new() -> Self {
633 Self {
634 v: -65.0,
635 h_na: 0.6,
636 n_k: 0.3,
637 m_t: 0.0,
638 h_t: 1.0,
639 g_na: 100.0,
640 g_k: 10.0,
641 g_t: 2.0,
642 g_l: 0.05,
643 e_na: 50.0,
644 e_k: -90.0,
645 e_ca: 120.0,
646 e_l: -70.0,
647 dt: 0.02,
648 v_threshold: -20.0,
649 }
650 }
651 pub fn step(&mut self, current: f64) -> i32 {
652 let v_prev = self.v;
653 for _ in 0..5 {
654 let m_na = 1.0 / (1.0 + (-(self.v + 37.0) / 7.0).exp());
655 let h_na_inf = 1.0 / (1.0 + ((self.v + 41.0) / 4.0).exp());
656 let n_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 12.0).exp());
657 let m_t_inf = 1.0 / (1.0 + (-(self.v + 57.0) / 6.5).exp());
658 let h_t_inf = 1.0 / (1.0 + ((self.v + 81.0) / 4.0).exp());
659 let tau_h_na = (1.0
661 / (0.128 * (-(self.v + 46.0) / 18.0).exp()
662 + 4.0 / (1.0 + (-(self.v + 23.0) / 5.0).exp())))
663 .max(0.1);
664 let tau_n_k = (1.0 / (0.032 * 5.0 + 0.5 * (-(self.v + 40.0) / 40.0).exp())).max(0.1);
665 let tau_h_t = if self.v < -81.0 {
666 (30.8
667 + 211.4 * ((self.v + 115.2) / 5.0).exp()
668 / (1.0 + ((self.v + 86.0) / 3.2).exp()))
669 .max(0.1)
670 } else {
671 10.0
672 };
673 self.h_na += (h_na_inf - self.h_na) / tau_h_na * self.dt;
674 self.n_k += (n_inf - self.n_k) / tau_n_k * self.dt;
675 self.m_t = m_t_inf; self.h_t += (h_t_inf - self.h_t) / tau_h_t * self.dt;
677 let i_na = self.g_na * m_na.powi(3) * self.h_na * (self.v - self.e_na);
678 let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
679 let i_t = self.g_t * self.m_t.powi(2) * self.h_t * (self.v - self.e_ca);
680 let i_l = self.g_l * (self.v - self.e_l);
681 self.v += (-i_na - i_k - i_t - i_l + current) * self.dt;
682 }
683 if self.v >= self.v_threshold && v_prev < self.v_threshold {
684 1
685 } else {
686 0
687 }
688 }
689 pub fn reset(&mut self) {
690 self.v = -65.0;
691 self.h_na = 0.6;
692 self.n_k = 0.3;
693 self.m_t = 0.0;
694 self.h_t = 1.0;
695 }
696}
697impl Default for DestexheThalamicNeuron {
698 fn default() -> Self {
699 Self::new()
700 }
701}
702
703#[derive(Clone, Debug)]
705pub struct HuberBraunNeuron {
706 pub v: f64,
707 pub a_sd: f64,
708 pub a_sr: f64,
709 pub g_sd: f64,
710 pub g_sr: f64,
711 pub g_l: f64,
712 pub e_sd: f64,
713 pub e_sr: f64,
714 pub e_l: f64,
715 pub tau_sd: f64,
716 pub tau_sr: f64,
717 pub eta: f64,
718 pub dt: f64,
719 pub v_threshold: f64,
720}
721
722impl HuberBraunNeuron {
723 pub fn new() -> Self {
724 Self {
725 v: -50.0,
726 a_sd: 0.0,
727 a_sr: 0.0,
728 g_sd: 1.5,
729 g_sr: 0.4,
730 g_l: 0.1,
731 e_sd: 50.0,
732 e_sr: -90.0,
733 e_l: -60.0,
734 tau_sd: 10.0,
735 tau_sr: 20.0,
736 eta: 0.012,
737 dt: 0.1,
738 v_threshold: -20.0,
739 }
740 }
741 pub fn step(&mut self, current: f64) -> i32 {
742 let v_prev = self.v;
743 let sd_inf = 1.0 / (1.0 + (-(self.v + 40.0) / 6.0).exp());
745 let sr_inf = 1.0 / (1.0 + ((self.v + 40.0) / 6.0).exp());
746 self.a_sd += (sd_inf - self.a_sd) / self.tau_sd * self.dt;
747 self.a_sr += (sr_inf - self.a_sr) / self.tau_sr * self.dt;
748 let i_sd = self.g_sd * self.a_sd * (self.v - self.e_sd);
749 let i_sr = self.g_sr * self.a_sr * (self.v - self.e_sr);
750 let i_l = self.g_l * (self.v - self.e_l);
751 self.v += (-i_sd - i_sr - i_l + current) * self.dt;
752 if self.v >= self.v_threshold && v_prev < self.v_threshold {
753 1
754 } else {
755 0
756 }
757 }
758 pub fn reset(&mut self) {
759 self.v = -50.0;
760 self.a_sd = 0.0;
761 self.a_sr = 0.0;
762 }
763}
764impl Default for HuberBraunNeuron {
765 fn default() -> Self {
766 Self::new()
767 }
768}
769
770#[derive(Clone, Debug)]
772pub struct GolombFSNeuron {
773 pub v: f64,
774 pub h: f64,
775 pub n: f64,
776 pub p: f64,
777 pub g_na: f64,
778 pub g_k: f64,
779 pub g_kv3: f64,
780 pub g_l: f64,
781 pub e_na: f64,
782 pub e_k: f64,
783 pub e_l: f64,
784 pub dt: f64,
785 pub v_threshold: f64,
786}
787
788impl GolombFSNeuron {
789 pub fn new() -> Self {
790 Self {
791 v: -65.0,
792 h: 0.9,
793 n: 0.1,
794 p: 0.0,
795 g_na: 112.5,
796 g_k: 225.0,
797 g_kv3: 150.0,
798 g_l: 0.25,
799 e_na: 50.0,
800 e_k: -90.0,
801 e_l: -70.0,
802 dt: 0.01,
803 v_threshold: -20.0,
804 }
805 }
806 fn derivatives(&self, v: f64, h: f64, n: f64, p: f64, current: f64) -> [f64; 4] {
811 let m_inf = 1.0 / (1.0 + (-(v + 24.0) / 11.5).exp());
812 let h_inf = 1.0 / (1.0 + ((v + 58.3) / 6.7).exp());
813 let n_inf = 1.0 / (1.0 + (-(v + 12.4) / 6.8).exp());
814 let p_inf = 1.0 / (1.0 + (-(v + 3.0) / 8.0).exp());
815 let tau_h = 0.5 + 14.0 / (1.0 + ((v + 60.0) / 12.0).exp());
816 let tau_n = 0.087 + 11.4 / (1.0 + ((v + 14.6) / 8.6).exp());
817 let tau_p = 0.1 + 4.0 / (1.0 + ((v + 25.0) / 10.0).exp());
818 let dh = (h_inf - h) / tau_h;
819 let dn = (n_inf - n) / tau_n;
820 let dp = (p_inf - p) / tau_p;
821 let i_na = self.g_na * m_inf * m_inf * m_inf * h * (v - self.e_na);
822 let i_k = self.g_k * n * n * n * n * (v - self.e_k);
823 let i_kv3 = self.g_kv3 * p * p * (v - self.e_k);
824 let i_l = self.g_l * (v - self.e_l);
825 let dv = -i_na - i_k - i_kv3 - i_l + current;
826 [dv, dh, dn, dp]
827 }
828
829 fn rk4_substep(&self, s: [f64; 4], current: f64) -> [f64; 4] {
832 let dt = self.dt;
833 let k1 = self.derivatives(s[0], s[1], s[2], s[3], current);
834 let k2 = self.derivatives(
835 s[0] + 0.5 * dt * k1[0],
836 s[1] + 0.5 * dt * k1[1],
837 s[2] + 0.5 * dt * k1[2],
838 s[3] + 0.5 * dt * k1[3],
839 current,
840 );
841 let k3 = self.derivatives(
842 s[0] + 0.5 * dt * k2[0],
843 s[1] + 0.5 * dt * k2[1],
844 s[2] + 0.5 * dt * k2[2],
845 s[3] + 0.5 * dt * k2[3],
846 current,
847 );
848 let k4 = self.derivatives(
849 s[0] + dt * k3[0],
850 s[1] + dt * k3[1],
851 s[2] + dt * k3[2],
852 s[3] + dt * k3[3],
853 current,
854 );
855 let mut out = [0.0_f64; 4];
856 for i in 0..4 {
857 out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
858 }
859 out
860 }
861
862 pub fn step(&mut self, current: f64) -> i32 {
863 let v_prev = self.v;
864 let mut s = [self.v, self.h, self.n, self.p];
865 for _ in 0..10 {
866 s = self.rk4_substep(s, current);
867 }
868 self.v = s[0];
869 self.h = s[1];
870 self.n = s[2];
871 self.p = s[3];
872 if self.v >= self.v_threshold && v_prev < self.v_threshold {
873 1
874 } else {
875 0
876 }
877 }
878 pub fn reset(&mut self) {
879 self.v = -65.0;
880 self.h = 0.9;
881 self.n = 0.1;
882 self.p = 0.0;
883 }
884}
885impl Default for GolombFSNeuron {
886 fn default() -> Self {
887 Self::new()
888 }
889}
890
891#[derive(Clone, Debug)]
893pub struct PospischilNeuron {
894 pub v: f64,
895 pub m: f64,
896 pub h: f64,
897 pub n: f64,
898 pub p: f64,
899 pub g_na: f64,
900 pub g_k: f64,
901 pub g_m: f64,
902 pub g_l: f64,
903 pub e_na: f64,
904 pub e_k: f64,
905 pub e_l: f64,
906 pub c_m: f64,
907 pub vt: f64,
908 pub dt: f64,
909 pub v_threshold: f64,
910}
911
912impl PospischilNeuron {
913 pub fn new() -> Self {
914 Self {
915 v: -70.0,
916 m: 0.05,
917 h: 0.6,
918 n: 0.3,
919 p: 0.0,
920 g_na: 50.0,
921 g_k: 5.0,
922 g_m: 0.07,
923 g_l: 0.1,
924 e_na: 50.0,
925 e_k: -90.0,
926 e_l: -70.0,
927 c_m: 1.0,
928 vt: -56.2,
929 dt: 0.025,
930 v_threshold: -20.0,
931 }
932 }
933 fn derivatives(&self, v: f64, m: f64, h: f64, n: f64, p: f64, current: f64) -> [f64; 5] {
938 let dv_vt = v - self.vt;
939 let x_m = dv_vt - 13.0;
940 let am = if x_m.abs() < 1e-6 {
941 1.28
942 } else {
943 -0.32 * x_m / ((-(x_m) / 4.0).exp() - 1.0)
944 };
945 let x_bm = dv_vt - 40.0;
946 let bm = if x_bm.abs() < 1e-6 {
947 1.4
948 } else {
949 0.28 * x_bm / ((x_bm / 5.0).exp() - 1.0)
950 };
951 let ah = 0.128 * (-(dv_vt - 17.0) / 18.0).exp();
952 let bh = 4.0 / (1.0 + (-(dv_vt - 40.0) / 5.0).exp());
953 let x_n = dv_vt - 15.0;
954 let an = if x_n.abs() < 1e-6 {
955 0.16
956 } else {
957 -0.032 * x_n / ((-(x_n) / 5.0).exp() - 1.0)
958 };
959 let bn = 0.5 * (-(dv_vt - 10.0) / 40.0).exp();
960 let p_inf = 1.0 / (1.0 + (-(v + 35.0) / 10.0).exp());
961 let tau_p = 608.0 / (3.3 * ((v + 35.0) / 20.0).exp() + (-(v + 35.0) / 20.0).exp());
962 let dm = am * (1.0 - m) - bm * m;
963 let dh = ah * (1.0 - h) - bh * h;
964 let dn = an * (1.0 - n) - bn * n;
965 let dp = (p_inf - p) / tau_p;
966 let i_na = self.g_na * m * m * m * h * (v - self.e_na);
967 let i_k = self.g_k * n * n * n * n * (v - self.e_k);
968 let i_m = self.g_m * p * (v - self.e_k);
969 let i_l = self.g_l * (v - self.e_l);
970 let dv = (-i_na - i_k - i_m - i_l + current) / self.c_m;
971 [dv, dm, dh, dn, dp]
972 }
973
974 fn rk4_substep(&self, s: [f64; 5], current: f64) -> [f64; 5] {
977 let dt = self.dt;
978 let k1 = self.derivatives(s[0], s[1], s[2], s[3], s[4], current);
979 let k2 = self.derivatives(
980 s[0] + 0.5 * dt * k1[0],
981 s[1] + 0.5 * dt * k1[1],
982 s[2] + 0.5 * dt * k1[2],
983 s[3] + 0.5 * dt * k1[3],
984 s[4] + 0.5 * dt * k1[4],
985 current,
986 );
987 let k3 = self.derivatives(
988 s[0] + 0.5 * dt * k2[0],
989 s[1] + 0.5 * dt * k2[1],
990 s[2] + 0.5 * dt * k2[2],
991 s[3] + 0.5 * dt * k2[3],
992 s[4] + 0.5 * dt * k2[4],
993 current,
994 );
995 let k4 = self.derivatives(
996 s[0] + dt * k3[0],
997 s[1] + dt * k3[1],
998 s[2] + dt * k3[2],
999 s[3] + dt * k3[3],
1000 s[4] + dt * k3[4],
1001 current,
1002 );
1003 let mut out = [0.0_f64; 5];
1004 for i in 0..5 {
1005 out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
1006 }
1007 out
1008 }
1009
1010 pub fn step(&mut self, current: f64) -> i32 {
1011 let v_prev = self.v;
1012 let mut s = [self.v, self.m, self.h, self.n, self.p];
1013 for _ in 0..4 {
1014 s = self.rk4_substep(s, current);
1015 }
1016 self.v = s[0];
1017 self.m = s[1];
1018 self.h = s[2];
1019 self.n = s[3];
1020 self.p = s[4];
1021 if self.v >= self.v_threshold && v_prev < self.v_threshold {
1022 1
1023 } else {
1024 0
1025 }
1026 }
1027 pub fn reset(&mut self) {
1028 self.v = -70.0;
1029 self.m = 0.05;
1030 self.h = 0.6;
1031 self.n = 0.3;
1032 self.p = 0.0;
1033 }
1034}
1035impl Default for PospischilNeuron {
1036 fn default() -> Self {
1037 Self::new()
1038 }
1039}
1040
1041#[derive(Clone, Debug)]
1043pub struct MainenSejnowskiNeuron {
1044 pub vs: f64,
1045 pub va: f64,
1046 pub m: f64,
1047 pub h: f64,
1048 pub n: f64,
1049 pub kappa: f64,
1050 pub g_na: f64,
1051 pub g_k: f64,
1052 pub g_l: f64,
1053 pub e_na: f64,
1054 pub e_k: f64,
1055 pub e_l: f64,
1056 pub c_s: f64,
1057 pub c_a: f64,
1058 pub dt: f64,
1059 pub v_threshold: f64,
1060}
1061
1062impl MainenSejnowskiNeuron {
1063 pub fn new() -> Self {
1064 Self {
1065 vs: -65.0,
1066 va: -65.0,
1067 m: 0.05,
1068 h: 0.6,
1069 n: 0.3,
1070 kappa: 10.0,
1071 g_na: 3000.0,
1072 g_k: 1500.0,
1073 g_l: 1.0,
1074 e_na: 50.0,
1075 e_k: -90.0,
1076 e_l: -70.0,
1077 c_s: 1.0,
1078 c_a: 0.1,
1079 dt: 0.005,
1080 v_threshold: -20.0,
1081 }
1082 }
1083 pub fn step(&mut self, current: f64) -> i32 {
1084 let v_prev = self.vs;
1085 for _ in 0..20 {
1086 let x_am = self.va + 25.0;
1088 let am = if x_am.abs() < 1e-6 {
1089 0.182 * 9.0
1090 } else {
1091 0.182 * x_am / (1.0 - (-(x_am) / 9.0).exp() + 1e-12)
1092 };
1093 let bm = if x_am.abs() < 1e-6 {
1094 0.124 * 9.0
1095 } else {
1096 -0.124 * x_am / (1.0 - ((x_am) / 9.0).exp() + 1e-12)
1097 };
1098 let x_ah = self.va + 40.0;
1099 let ah = if x_ah.abs() < 1e-6 {
1100 0.024 * 5.0
1101 } else {
1102 0.024 * x_ah / (1.0 - (-(x_ah) / 5.0).exp() + 1e-12)
1103 };
1104 let x_bh = self.va + 65.0;
1105 let bh = if x_bh.abs() < 1e-6 {
1106 0.0091 * 5.0
1107 } else {
1108 -0.0091 * x_bh / (1.0 - ((x_bh) / 5.0).exp() + 1e-12)
1109 };
1110 let x_an = self.va - 20.0;
1111 let an = if x_an.abs() < 1e-6 {
1112 0.02 * 9.0
1113 } else {
1114 0.02 * x_an / (1.0 - (-(x_an) / 9.0).exp() + 1e-12)
1115 };
1116 let bn = if x_an.abs() < 1e-6 {
1117 0.002 * 9.0
1118 } else {
1119 -0.002 * x_an / (1.0 - ((x_an) / 9.0).exp() + 1e-12)
1120 };
1121 self.m = (self.m + (am * (1.0 - self.m) - bm * self.m) * self.dt).clamp(0.0, 1.0);
1122 self.h = (self.h + (ah * (1.0 - self.h) - bh * self.h) * self.dt).clamp(0.0, 1.0);
1123 self.n = (self.n + (an * (1.0 - self.n) - bn * self.n) * self.dt).clamp(0.0, 1.0);
1124 let i_na = self.g_na * self.m.powi(3) * self.h * (self.va - self.e_na);
1125 let i_k = self.g_k * self.n * (self.va - self.e_k);
1126 let i_l_s = self.g_l * (self.vs - self.e_l);
1127 self.vs = (self.vs
1128 + (-i_l_s + self.kappa * (self.va - self.vs) + current) / self.c_s * self.dt)
1129 .clamp(-200.0, 200.0);
1130 self.va = (self.va
1131 + (-i_na - i_k + self.kappa * (self.vs - self.va)) / self.c_a * self.dt)
1132 .clamp(-200.0, 200.0);
1133 }
1134 if self.vs >= self.v_threshold && v_prev < self.v_threshold {
1135 1
1136 } else {
1137 0
1138 }
1139 }
1140 pub fn reset(&mut self) {
1141 self.vs = -65.0;
1142 self.va = -65.0;
1143 self.m = 0.05;
1144 self.h = 0.6;
1145 self.n = 0.3;
1146 }
1147}
1148impl Default for MainenSejnowskiNeuron {
1149 fn default() -> Self {
1150 Self::new()
1151 }
1152}
1153
1154#[derive(Clone, Debug)]
1159pub struct DeSchutterPurkinjeNeuron {
1160 pub v: f64,
1161 pub h_na: f64,
1162 pub n_k: f64,
1163 pub m_cap: f64,
1164 pub h_cap: f64,
1165 pub q_kca: f64,
1166 pub ca: f64,
1167 pub g_na: f64,
1168 pub g_k: f64,
1169 pub g_cap: f64,
1170 pub g_kca: f64,
1171 pub g_l: f64,
1172 pub e_na: f64,
1173 pub e_k: f64,
1174 pub e_ca: f64,
1175 pub e_l: f64,
1176 pub ca_decay: f64,
1177 pub f_ca: f64,
1178 pub dt: f64,
1179 pub v_threshold: f64,
1180}
1181
1182impl DeSchutterPurkinjeNeuron {
1183 pub fn new() -> Self {
1184 Self {
1185 v: -68.0,
1186 h_na: 0.8,
1187 n_k: 0.1,
1188 m_cap: 0.0,
1189 h_cap: 0.9,
1190 q_kca: 0.0,
1191 ca: 0.0001,
1192 g_na: 125.0,
1193 g_k: 10.0,
1194 g_cap: 45.0,
1195 g_kca: 35.0,
1196 g_l: 0.5,
1197 e_na: 45.0,
1198 e_k: -85.0,
1199 e_ca: 135.0,
1200 e_l: -68.0,
1201 ca_decay: 0.02,
1202 f_ca: 0.00024,
1203 dt: 0.01,
1204 v_threshold: -20.0,
1205 }
1206 }
1207 fn valid(&self) -> bool {
1208 self.v.is_finite()
1209 && self.h_na.is_finite()
1210 && self.n_k.is_finite()
1211 && self.m_cap.is_finite()
1212 && self.h_cap.is_finite()
1213 && self.q_kca.is_finite()
1214 && self.ca.is_finite()
1215 && self.ca >= 0.0
1216 && self.g_na.is_finite()
1217 && self.g_na >= 0.0
1218 && self.g_k.is_finite()
1219 && self.g_k >= 0.0
1220 && self.g_cap.is_finite()
1221 && self.g_cap >= 0.0
1222 && self.g_kca.is_finite()
1223 && self.g_kca >= 0.0
1224 && self.g_l.is_finite()
1225 && self.g_l >= 0.0
1226 && self.e_na.is_finite()
1227 && self.e_k.is_finite()
1228 && self.e_ca.is_finite()
1229 && self.e_l.is_finite()
1230 && self.ca_decay.is_finite()
1231 && self.ca_decay >= 0.0
1232 && self.f_ca.is_finite()
1233 && self.f_ca >= 0.0
1234 && self.dt.is_finite()
1235 && self.dt > 0.0
1236 && self.v_threshold.is_finite()
1237 }
1238 fn derivatives(&self, s: [f64; 7], current: f64) -> [f64; 7] {
1239 let v = s[0];
1240 let h_na = s[1];
1241 let n_k = s[2];
1242 let m_cap = s[3];
1243 let h_cap = s[4];
1244 let q_kca = s[5];
1245 let ca = s[6].max(0.0);
1246 let m_na = 1.0 / (1.0 + (-(v + 35.0) / 7.5).exp());
1247 let h_na_inf = 1.0 / (1.0 + ((v + 55.0) / 7.0).exp());
1248 let n_inf = 1.0 / (1.0 + (-(v + 30.0) / 15.0).exp());
1249 let m_cap_inf = 1.0 / (1.0 + (-(v + 19.0) / 5.5).exp());
1250 let h_cap_inf = 1.0 / (1.0 + ((v + 48.0) / 7.0).exp());
1251 let q_inf = ca / (ca + 0.0002);
1252 let tau_h_na = 0.5 + 14.0 / (1.0 + ((v + 40.0) / 12.0).exp());
1253 let tau_n_k = 1.0 + 11.0 / (1.0 + ((v + 15.0) / 8.0).exp());
1254 let i_na = self.g_na * m_na * m_na * m_na * h_na * (v - self.e_na);
1255 let i_k = self.g_k * n_k * n_k * n_k * n_k * (v - self.e_k);
1256 let i_cap = self.g_cap * m_cap * m_cap * h_cap * (v - self.e_ca);
1257 let i_kca = self.g_kca * q_kca * (v - self.e_k);
1258 let i_l = self.g_l * (v - self.e_l);
1259 [
1260 -i_na - i_k - i_cap - i_kca - i_l + current,
1261 (h_na_inf - h_na) / tau_h_na,
1262 (n_inf - n_k) / tau_n_k,
1263 (m_cap_inf - m_cap) / 0.3,
1264 (h_cap_inf - h_cap) / 45.0,
1265 q_inf - q_kca,
1266 -self.f_ca * i_cap - self.ca_decay * ca,
1267 ]
1268 }
1269 fn rk4_substep(&self, s: [f64; 7], current: f64) -> [f64; 7] {
1270 let dt = self.dt;
1271 let k1 = self.derivatives(s, current);
1272 let mut s2 = [0.0; 7];
1273 let mut s3 = [0.0; 7];
1274 let mut s4 = [0.0; 7];
1275 for i in 0..7 {
1276 s2[i] = s[i] + 0.5 * dt * k1[i];
1277 }
1278 let k2 = self.derivatives(s2, current);
1279 for i in 0..7 {
1280 s3[i] = s[i] + 0.5 * dt * k2[i];
1281 }
1282 let k3 = self.derivatives(s3, current);
1283 for i in 0..7 {
1284 s4[i] = s[i] + dt * k3[i];
1285 }
1286 let k4 = self.derivatives(s4, current);
1287 let mut next = [0.0; 7];
1288 for i in 0..7 {
1289 next[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
1290 }
1291 next[6] = next[6].max(0.0);
1292 next
1293 }
1294 pub fn step(&mut self, current: f64) -> i32 {
1295 if !current.is_finite() || !self.valid() {
1296 return 0;
1297 }
1298 let v_prev = self.v;
1299 let mut state = [
1300 self.v, self.h_na, self.n_k, self.m_cap, self.h_cap, self.q_kca, self.ca,
1301 ];
1302 for _ in 0..5 {
1303 state = self.rk4_substep(state, current);
1304 if !state.iter().all(|value| value.is_finite()) {
1305 return 0;
1306 }
1307 }
1308 self.v = state[0];
1309 self.h_na = state[1];
1310 self.n_k = state[2];
1311 self.m_cap = state[3];
1312 self.h_cap = state[4];
1313 self.q_kca = state[5];
1314 self.ca = state[6];
1315 if self.v >= self.v_threshold && v_prev < self.v_threshold {
1316 1
1317 } else {
1318 0
1319 }
1320 }
1321 pub fn reset(&mut self) {
1322 self.v = -68.0;
1323 self.h_na = 0.8;
1324 self.n_k = 0.1;
1325 self.m_cap = 0.0;
1326 self.h_cap = 0.9;
1327 self.q_kca = 0.0;
1328 self.ca = 0.0001;
1329 }
1330}
1331impl Default for DeSchutterPurkinjeNeuron {
1332 fn default() -> Self {
1333 Self::new()
1334 }
1335}
1336
1337#[derive(Clone, Debug)]
1339pub struct PlantR15Neuron {
1340 pub v: f64,
1341 pub m: f64,
1342 pub h: f64,
1343 pub n: f64,
1344 pub ca: f64,
1345 pub g_na: f64,
1346 pub g_k: f64,
1347 pub g_ca: f64,
1348 pub g_l: f64,
1349 pub g_kca: f64,
1350 pub e_na: f64,
1351 pub e_k: f64,
1352 pub e_ca: f64,
1353 pub e_l: f64,
1354 pub c_m: f64,
1355 pub k_ca: f64,
1356 pub tau_ca: f64,
1357 pub dt: f64,
1358 pub v_threshold: f64,
1359}
1360
1361impl PlantR15Neuron {
1362 pub fn new() -> Self {
1363 Self {
1364 v: -50.0,
1365 m: 0.05,
1366 h: 0.6,
1367 n: 0.3,
1368 ca: 0.1,
1369 g_na: 4.0,
1370 g_k: 0.3,
1371 g_ca: 0.004,
1372 g_l: 0.003,
1373 g_kca: 0.03,
1374 e_na: 30.0,
1375 e_k: -75.0,
1376 e_ca: 140.0,
1377 e_l: -40.0,
1378 c_m: 1.0,
1379 k_ca: 0.0085,
1380 tau_ca: 500.0,
1381 dt: 0.05,
1382 v_threshold: -10.0,
1383 }
1384 }
1385 pub fn step(&mut self, current: f64) -> i32 {
1386 let v_prev = self.v;
1387 for _ in 0..5 {
1388 let am = safe_rate(0.1, 50.0, self.v, 10.0, 1.0);
1389 let bm = 4.0 * (-(self.v + 75.0) / 18.0).exp();
1390 let ah = 0.07 * (-(self.v + 50.0) / 20.0).exp();
1391 let bh = 1.0 / (1.0 + (-(self.v + 20.0) / 10.0).exp());
1392 let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
1393 let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
1394 self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
1395 self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
1396 self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
1397 let m_ca = 1.0 / (1.0 + (-(self.v + 25.0) / 5.0).exp());
1398 let kca_act = self.ca / (0.5 + self.ca);
1399 let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
1400 let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
1401 let i_ca = self.g_ca * m_ca.powi(2) * (self.v - self.e_ca);
1402 let i_kca = self.g_kca * kca_act * (self.v - self.e_k);
1403 let i_l = self.g_l * (self.v - self.e_l);
1404 self.v += (-i_na - i_k - i_ca - i_kca - i_l + current) / self.c_m * self.dt;
1405 self.ca = (self.ca + (-self.k_ca * i_ca - self.ca / self.tau_ca) * self.dt).max(0.0);
1406 }
1407 if self.v >= self.v_threshold && v_prev < self.v_threshold {
1408 1
1409 } else {
1410 0
1411 }
1412 }
1413 pub fn reset(&mut self) {
1414 self.v = -50.0;
1415 self.m = 0.05;
1416 self.h = 0.6;
1417 self.n = 0.3;
1418 self.ca = 0.1;
1419 }
1420}
1421impl Default for PlantR15Neuron {
1422 fn default() -> Self {
1423 Self::new()
1424 }
1425}
1426
1427#[derive(Clone, Debug)]
1429pub struct PrescottNeuron {
1430 pub v: f64,
1431 pub w: f64,
1432 pub g_fast: f64,
1433 pub g_slow: f64,
1434 pub g_l: f64,
1435 pub e_fast: f64,
1436 pub e_slow: f64,
1437 pub e_l: f64,
1438 pub beta_w: f64,
1439 pub gamma_w: f64,
1440 pub tau_w: f64,
1441 pub phi: f64,
1442 pub dt: f64,
1443 pub v_threshold: f64,
1444}
1445
1446impl PrescottNeuron {
1447 pub fn new() -> Self {
1448 Self {
1449 v: -65.0,
1450 w: 0.0,
1451 g_fast: 20.0,
1452 g_slow: 20.0,
1453 g_l: 2.0,
1454 e_fast: 50.0,
1455 e_slow: -100.0,
1456 e_l: -70.0,
1457 beta_w: -21.0,
1458 gamma_w: 15.0,
1459 tau_w: 100.0,
1460 phi: 0.15,
1461 dt: 0.1,
1462 v_threshold: -20.0,
1463 }
1464 }
1465
1466 fn sigmoid(x: f64) -> f64 {
1467 if x >= 0.0 {
1468 let z = (-x).exp();
1469 1.0 / (1.0 + z)
1470 } else {
1471 let z = x.exp();
1472 z / (1.0 + z)
1473 }
1474 }
1475
1476 fn valid_state(v: f64, w: f64) -> bool {
1477 v.is_finite() && w.is_finite() && (0.0..=1.0).contains(&w)
1478 }
1479
1480 fn valid_runtime(&self) -> bool {
1481 Self::valid_state(self.v, self.w)
1482 && self.g_fast.is_finite()
1483 && self.g_fast >= 0.0
1484 && self.g_slow.is_finite()
1485 && self.g_slow >= 0.0
1486 && self.g_l.is_finite()
1487 && self.g_l >= 0.0
1488 && self.e_fast.is_finite()
1489 && self.e_slow.is_finite()
1490 && self.e_l.is_finite()
1491 && self.beta_w.is_finite()
1492 && self.gamma_w.is_finite()
1493 && self.gamma_w > 0.0
1494 && self.tau_w.is_finite()
1495 && self.tau_w > 0.0
1496 && self.phi.is_finite()
1497 && self.phi >= 0.0
1498 && self.dt.is_finite()
1499 && self.dt > 0.0
1500 && self.v_threshold.is_finite()
1501 }
1502
1503 fn derivatives(&self, v: f64, w: f64, current: f64) -> Option<(f64, f64)> {
1504 if !Self::valid_state(v, w) {
1505 return None;
1506 }
1507 let m_inf = Self::sigmoid((v + 20.0) / 15.0);
1508 let w_inf = Self::sigmoid((v - self.beta_w) / self.gamma_w);
1509 let i_fast = self.g_fast * m_inf * (v - self.e_fast);
1510 let i_slow = self.g_slow * w * (v - self.e_slow);
1511 let i_l = self.g_l * (v - self.e_l);
1512 let dv = -i_fast - i_slow - i_l + current;
1513 let dw = self.phi * (w_inf - w) / self.tau_w;
1514 if dv.is_finite() && dw.is_finite() {
1515 Some((dv, dw))
1516 } else {
1517 None
1518 }
1519 }
1520
1521 fn rk4_step(&self, current: f64) -> Option<(f64, f64)> {
1522 let dt = self.dt;
1523 let (k1_v, k1_w) = self.derivatives(self.v, self.w, current)?;
1524 let (k2_v, k2_w) =
1525 self.derivatives(self.v + 0.5 * dt * k1_v, self.w + 0.5 * dt * k1_w, current)?;
1526 let (k3_v, k3_w) =
1527 self.derivatives(self.v + 0.5 * dt * k2_v, self.w + 0.5 * dt * k2_w, current)?;
1528 let (k4_v, k4_w) = self.derivatives(self.v + dt * k3_v, self.w + dt * k3_w, current)?;
1529 let next_v = self.v + dt * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v) / 6.0;
1530 let next_w = self.w + dt * (k1_w + 2.0 * k2_w + 2.0 * k3_w + k4_w) / 6.0;
1531 if Self::valid_state(next_v, next_w) {
1532 Some((next_v, next_w))
1533 } else {
1534 None
1535 }
1536 }
1537
1538 pub fn step(&mut self, current: f64) -> i32 {
1539 if !current.is_finite() || !self.valid_runtime() {
1540 return 0;
1541 }
1542 let v_prev = self.v;
1543 let Some((next_v, next_w)) = self.rk4_step(current) else {
1544 return 0;
1545 };
1546 self.v = next_v;
1547 self.w = next_w;
1548 if self.v >= self.v_threshold && v_prev < self.v_threshold {
1549 1
1550 } else {
1551 0
1552 }
1553 }
1554
1555 pub fn reset(&mut self) {
1556 self.v = -65.0;
1557 self.w = 0.0;
1558 }
1559}
1560impl Default for PrescottNeuron {
1561 fn default() -> Self {
1562 Self::new()
1563 }
1564}
1565
1566#[derive(Clone, Debug)]
1568pub struct MihalasNieburNeuron {
1569 pub v: f64,
1570 pub theta: f64,
1571 pub i1: f64,
1572 pub i2: f64,
1573 pub v_rest: f64,
1574 pub v_reset: f64,
1575 pub theta_reset: f64,
1576 pub theta_inf: f64,
1577 pub tau_v: f64,
1578 pub tau_theta: f64,
1579 pub tau_1: f64,
1580 pub tau_2: f64,
1581 pub a: f64,
1582 pub b: f64,
1583 pub r1: f64,
1584 pub r2: f64,
1585 pub dt: f64,
1586}
1587
1588impl MihalasNieburNeuron {
1589 pub fn new() -> Self {
1590 Self {
1591 v: 0.0,
1592 theta: 1.0,
1593 i1: 0.0,
1594 i2: 0.0,
1595 v_rest: 0.0,
1596 v_reset: 0.0,
1597 theta_reset: 1.0,
1598 theta_inf: 1.0,
1599 tau_v: 10.0,
1600 tau_theta: 100.0,
1601 tau_1: 10.0,
1602 tau_2: 200.0,
1603 a: 0.0,
1604 b: 0.0,
1605 r1: 0.0,
1606 r2: 0.0,
1607 dt: 1.0,
1608 }
1609 }
1610
1611 fn finite_values(values: &[f64]) -> bool {
1612 values.iter().all(|value| value.is_finite())
1613 }
1614
1615 fn valid_runtime(&self) -> bool {
1616 Self::finite_values(&[
1617 self.v,
1618 self.theta,
1619 self.i1,
1620 self.i2,
1621 self.v_rest,
1622 self.v_reset,
1623 self.theta_reset,
1624 self.theta_inf,
1625 self.tau_v,
1626 self.tau_theta,
1627 self.tau_1,
1628 self.tau_2,
1629 self.a,
1630 self.b,
1631 self.r1,
1632 self.r2,
1633 self.dt,
1634 ]) && self.tau_v > 0.0
1635 && self.tau_theta > 0.0
1636 && self.tau_1 > 0.0
1637 && self.tau_2 > 0.0
1638 && self.dt > 0.0
1639 }
1640
1641 fn derivatives(&self, v: f64, theta: f64, i1: f64, i2: f64, current: f64) -> [f64; 4] {
1642 [
1643 (-(v - self.v_rest) + i1 + i2 + current) / self.tau_v,
1644 (self.theta_inf - theta + self.a * (v - self.v_rest)) / self.tau_theta,
1645 -i1 / self.tau_1,
1646 -i2 / self.tau_2,
1647 ]
1648 }
1649
1650 fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
1651 [
1652 state[0] + scale * slope[0],
1653 state[1] + scale * slope[1],
1654 state[2] + scale * slope[2],
1655 state[3] + scale * slope[3],
1656 ]
1657 }
1658
1659 fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
1660 let state = [self.v, self.theta, self.i1, self.i2];
1661 let half_dt = 0.5 * self.dt;
1662 let k1 = self.derivatives(state[0], state[1], state[2], state[3], current);
1663 let s2 = Self::add_scaled(state, k1, half_dt);
1664 let k2 = self.derivatives(s2[0], s2[1], s2[2], s2[3], current);
1665 let s3 = Self::add_scaled(state, k2, half_dt);
1666 let k3 = self.derivatives(s3[0], s3[1], s3[2], s3[3], current);
1667 let s4 = Self::add_scaled(state, k3, self.dt);
1668 let k4 = self.derivatives(s4[0], s4[1], s4[2], s4[3], current);
1669 let candidate = [
1670 state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
1671 state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
1672 state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
1673 state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
1674 ];
1675 if Self::finite_values(&candidate) {
1676 Some(candidate)
1677 } else {
1678 None
1679 }
1680 }
1681
1682 pub fn step(&mut self, current: f64) -> i32 {
1683 if !current.is_finite() || !self.valid_runtime() {
1684 return 0;
1685 }
1686 let Some(candidate) = self.rk4_candidate(current) else {
1687 return 0;
1688 };
1689 self.v = candidate[0];
1690 self.theta = candidate[1];
1691 self.i1 = candidate[2];
1692 self.i2 = candidate[3];
1693 if self.v >= self.theta {
1694 self.v = self.v_reset + self.b * (self.v - self.v_rest);
1695 self.theta = self.theta_reset.max(self.theta);
1696 self.i1 += self.r1;
1697 self.i2 += self.r2;
1698 1
1699 } else {
1700 0
1701 }
1702 }
1703
1704 pub fn simulate(&mut self, n_steps: usize, current: f64) -> (Vec<f64>, i64) {
1710 let mut trace = Vec::with_capacity(n_steps);
1711 let mut spikes: i64 = 0;
1712 for _ in 0..n_steps {
1713 spikes += i64::from(self.step(current));
1714 trace.push(self.v);
1715 }
1716 (trace, spikes)
1717 }
1718
1719 pub fn reset(&mut self) {
1720 self.v = self.v_rest;
1721 self.theta = self.theta_reset;
1722 self.i1 = 0.0;
1723 self.i2 = 0.0;
1724 }
1725}
1726impl Default for MihalasNieburNeuron {
1727 fn default() -> Self {
1728 Self::new()
1729 }
1730}
1731
1732#[derive(Clone, Debug)]
1734pub struct GLIFNeuron {
1735 pub v: f64,
1736 pub theta: f64,
1737 pub theta_inf: f64,
1738 pub i_asc1: f64,
1739 pub i_asc2: f64,
1740 pub v_rest: f64,
1741 pub v_reset: f64,
1742 pub tau_m: f64,
1743 pub tau_theta: f64,
1744 pub tau_asc1: f64,
1745 pub tau_asc2: f64,
1746 pub a_theta: f64,
1747 pub delta_theta: f64,
1748 pub r_asc1: f64,
1749 pub r_asc2: f64,
1750 pub resistance: f64,
1751 pub dt: f64,
1752}
1753
1754impl GLIFNeuron {
1755 pub fn new() -> Self {
1756 Self {
1757 v: -70.0,
1758 theta: -50.0,
1759 theta_inf: -50.0,
1760 i_asc1: 0.0,
1761 i_asc2: 0.0,
1762 v_rest: -70.0,
1763 v_reset: -70.0,
1764 tau_m: 10.0,
1765 tau_theta: 100.0,
1766 tau_asc1: 10.0,
1767 tau_asc2: 200.0,
1768 a_theta: 0.01,
1769 delta_theta: 2.0,
1770 r_asc1: 1.0,
1771 r_asc2: 0.5,
1772 resistance: 1.0,
1773 dt: 1.0,
1774 }
1775 }
1776 fn finite_values(values: &[f64]) -> bool {
1777 values.iter().all(|value| value.is_finite())
1778 }
1779
1780 fn valid_runtime(&self) -> bool {
1781 Self::finite_values(&[
1782 self.v,
1783 self.theta,
1784 self.theta_inf,
1785 self.i_asc1,
1786 self.i_asc2,
1787 self.v_rest,
1788 self.v_reset,
1789 self.tau_m,
1790 self.tau_theta,
1791 self.tau_asc1,
1792 self.tau_asc2,
1793 self.a_theta,
1794 self.delta_theta,
1795 self.r_asc1,
1796 self.r_asc2,
1797 self.resistance,
1798 self.dt,
1799 ]) && self.tau_m > 0.0
1800 && self.tau_theta > 0.0
1801 && self.tau_asc1 > 0.0
1802 && self.tau_asc2 > 0.0
1803 && self.dt > 0.0
1804 && self.delta_theta >= 0.0
1805 && self.resistance >= 0.0
1806 }
1807
1808 fn derivatives(&self, v: f64, theta: f64, i_asc1: f64, i_asc2: f64, current: f64) -> [f64; 4] {
1809 [
1810 (-(v - self.v_rest) + self.resistance * current + i_asc1 + i_asc2) / self.tau_m,
1811 (self.theta_inf - theta + self.a_theta * (v - self.v_rest)) / self.tau_theta,
1812 -i_asc1 / self.tau_asc1,
1813 -i_asc2 / self.tau_asc2,
1814 ]
1815 }
1816
1817 fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
1818 [
1819 state[0] + scale * slope[0],
1820 state[1] + scale * slope[1],
1821 state[2] + scale * slope[2],
1822 state[3] + scale * slope[3],
1823 ]
1824 }
1825
1826 fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
1827 let state = [self.v, self.theta, self.i_asc1, self.i_asc2];
1828 let half_dt = 0.5 * self.dt;
1829 let k1 = self.derivatives(state[0], state[1], state[2], state[3], current);
1830 let s2 = Self::add_scaled(state, k1, half_dt);
1831 let k2 = self.derivatives(s2[0], s2[1], s2[2], s2[3], current);
1832 let s3 = Self::add_scaled(state, k2, half_dt);
1833 let k3 = self.derivatives(s3[0], s3[1], s3[2], s3[3], current);
1834 let s4 = Self::add_scaled(state, k3, self.dt);
1835 let k4 = self.derivatives(s4[0], s4[1], s4[2], s4[3], current);
1836 let candidate = [
1837 state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
1838 state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
1839 state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
1840 state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
1841 ];
1842 if Self::finite_values(&candidate) {
1843 Some(candidate)
1844 } else {
1845 None
1846 }
1847 }
1848
1849 pub fn step(&mut self, current: f64) -> i32 {
1850 if !current.is_finite() || !self.valid_runtime() {
1851 return 0;
1852 }
1853 let Some(candidate) = self.rk4_candidate(current) else {
1854 return 0;
1855 };
1856 self.v = candidate[0];
1857 self.theta = candidate[1];
1858 self.i_asc1 = candidate[2];
1859 self.i_asc2 = candidate[3];
1860 if self.v >= self.theta {
1861 self.v = self.v_reset;
1862 self.theta += self.delta_theta;
1863 self.i_asc1 += self.r_asc1;
1864 self.i_asc2 += self.r_asc2;
1865 1
1866 } else {
1867 0
1868 }
1869 }
1870
1871 pub fn simulate(&mut self, n_steps: usize, current: f64) -> (Vec<f64>, i64) {
1877 let mut trace = Vec::with_capacity(n_steps);
1878 let mut spikes: i64 = 0;
1879 for _ in 0..n_steps {
1880 spikes += i64::from(self.step(current));
1881 trace.push(self.v);
1882 }
1883 (trace, spikes)
1884 }
1885
1886 pub fn reset(&mut self) {
1887 self.v = self.v_rest;
1888 self.theta = self.theta_inf;
1889 self.i_asc1 = 0.0;
1890 self.i_asc2 = 0.0;
1891 }
1892}
1893impl Default for GLIFNeuron {
1894 fn default() -> Self {
1895 Self::new()
1896 }
1897}
1898
1899#[derive(Clone, Debug)]
1901pub struct GIFPopulationNeuron {
1902 pub v: f64,
1903 pub theta: f64,
1904 pub eta: f64,
1905 pub tau_m: f64,
1906 pub tau_eta: f64,
1907 pub delta_v: f64,
1908 pub lambda_0: f64,
1909 pub eta_increment: f64,
1910 pub v_rest: f64,
1911 pub v_reset: f64,
1912 pub dt: f64,
1913 pub seed: u64,
1914 rng: Xoshiro256PlusPlus,
1915}
1916
1917impl GIFPopulationNeuron {
1918 pub fn new(seed: u64) -> Self {
1919 Self {
1920 v: -65.0,
1921 theta: -50.0,
1922 eta: 0.0,
1923 tau_m: 20.0,
1924 tau_eta: 100.0,
1925 delta_v: 2.0,
1926 lambda_0: 0.001,
1927 eta_increment: 5.0,
1928 v_rest: -65.0,
1929 v_reset: -65.0,
1930 dt: 0.5,
1931 seed,
1932 rng: Xoshiro256PlusPlus::seed_from_u64(seed),
1933 }
1934 }
1935
1936 fn finite_values(values: &[f64]) -> bool {
1937 values.iter().all(|value| value.is_finite())
1938 }
1939
1940 fn valid_runtime(&self) -> bool {
1941 Self::finite_values(&[
1942 self.v,
1943 self.theta,
1944 self.eta,
1945 self.tau_m,
1946 self.tau_eta,
1947 self.delta_v,
1948 self.lambda_0,
1949 self.eta_increment,
1950 self.v_rest,
1951 self.v_reset,
1952 self.dt,
1953 ]) && self.tau_m > 0.0
1954 && self.tau_eta > 0.0
1955 && self.delta_v > 0.0
1956 && self.lambda_0 >= 0.0
1957 && self.dt > 0.0
1958 }
1959
1960 fn advance_subthreshold(&self, current: f64) -> Option<(f64, f64)> {
1961 let eta_decay = (-self.dt / self.tau_eta).exp();
1962 let membrane_decay = (-self.dt / self.tau_m).exp();
1963 let x0 = self.v - self.v_rest - current;
1964 let eta_new = self.eta * eta_decay;
1965 let x_new = if (self.tau_m - self.tau_eta).abs() <= 1e-12 {
1966 membrane_decay * (x0 - self.eta * self.dt / self.tau_m)
1967 } else {
1968 let coupling = self.tau_eta / (self.tau_eta - self.tau_m);
1969 x0 * membrane_decay - self.eta * coupling * (eta_decay - membrane_decay)
1970 };
1971 let v_new = self.v_rest + current + x_new;
1972 if Self::finite_values(&[v_new, eta_new]) {
1973 Some((v_new, eta_new))
1974 } else {
1975 None
1976 }
1977 }
1978
1979 fn spike_probability(&self, voltage: f64) -> f64 {
1980 if self.lambda_0 == 0.0 {
1981 return 0.0;
1982 }
1983 let exponent = ((voltage - self.theta) / self.delta_v).clamp(-745.0, 20.0);
1984 let hazard = self.lambda_0 * exponent.exp();
1985 (1.0 - (-hazard * self.dt).exp()).clamp(0.0, 1.0)
1986 }
1987
1988 pub fn step(&mut self, current: f64) -> i32 {
1989 if !current.is_finite() || !self.valid_runtime() {
1990 return 0;
1991 }
1992 let Some((v_candidate, eta_candidate)) = self.advance_subthreshold(current) else {
1993 return 0;
1994 };
1995 self.v = v_candidate;
1996 self.eta = eta_candidate;
1997 if self.rng.random::<f64>() < self.spike_probability(self.v) {
1998 self.v = self.v_reset;
1999 self.eta += self.eta_increment;
2000 1
2001 } else {
2002 0
2003 }
2004 }
2005
2006 pub fn reset(&mut self) {
2007 self.v = self.v_rest;
2008 self.eta = 0.0;
2009 self.rng = Xoshiro256PlusPlus::seed_from_u64(self.seed);
2010 }
2011}
2012
2013#[derive(Clone, Debug)]
2015pub struct AvRonCardiacNeuron {
2016 pub v: f64,
2017 pub h: f64,
2018 pub n: f64,
2019 pub s: f64,
2020 pub g_na: f64,
2021 pub g_k: f64,
2022 pub g_s: f64,
2023 pub g_l: f64,
2024 pub e_na: f64,
2025 pub e_k: f64,
2026 pub e_s: f64,
2027 pub e_l: f64,
2028 pub dt: f64,
2029 pub v_threshold: f64,
2030}
2031
2032impl AvRonCardiacNeuron {
2033 pub fn new() -> Self {
2034 Self {
2035 v: -60.0,
2036 h: 0.6,
2037 n: 0.3,
2038 s: 0.5,
2039 g_na: 80.0,
2040 g_k: 40.0,
2041 g_s: 20.0,
2042 g_l: 0.1,
2043 e_na: 40.0,
2044 e_k: -80.0,
2045 e_s: -25.0,
2046 e_l: -60.0,
2047 dt: 0.02,
2048 v_threshold: -20.0,
2049 }
2050 }
2051
2052 fn finite_values(values: &[f64]) -> bool {
2053 values.iter().all(|value| value.is_finite())
2054 }
2055
2056 fn gate_in_range(value: f64) -> bool {
2057 (0.0..=1.0).contains(&value)
2058 }
2059
2060 fn bounded_exp(value: f64) -> f64 {
2061 value.clamp(-745.0, 709.0).exp()
2062 }
2063
2064 fn sigmoid_pos(value: f64) -> f64 {
2065 1.0 / (1.0 + Self::bounded_exp(-value))
2066 }
2067
2068 fn sigmoid_neg(value: f64) -> f64 {
2069 1.0 / (1.0 + Self::bounded_exp(value))
2070 }
2071
2072 fn valid_runtime(&self) -> bool {
2073 Self::finite_values(&[
2074 self.v,
2075 self.h,
2076 self.n,
2077 self.s,
2078 self.g_na,
2079 self.g_k,
2080 self.g_s,
2081 self.g_l,
2082 self.e_na,
2083 self.e_k,
2084 self.e_s,
2085 self.e_l,
2086 self.dt,
2087 self.v_threshold,
2088 ]) && self.dt > 0.0
2089 && self.g_na >= 0.0
2090 && self.g_k >= 0.0
2091 && self.g_s >= 0.0
2092 && self.g_l >= 0.0
2093 && Self::gate_in_range(self.h)
2094 && Self::gate_in_range(self.n)
2095 && Self::gate_in_range(self.s)
2096 }
2097
2098 fn rates(&self, voltage: f64) -> [f64; 7] {
2099 [
2100 Self::sigmoid_pos((voltage + 40.0) / 7.0),
2101 Self::sigmoid_neg((voltage + 45.0) / 5.0),
2102 Self::sigmoid_pos((voltage + 40.0) / 15.0),
2103 Self::sigmoid_neg((voltage + 35.0) / 3.0),
2104 1.0 + 12.0 * Self::sigmoid_neg((voltage + 50.0) / 8.0),
2105 1.0 + 8.0 * Self::sigmoid_neg((voltage + 35.0) / 8.0),
2106 200.0 + 1000.0 * Self::sigmoid_neg((voltage + 30.0) / 5.0),
2107 ]
2108 }
2109
2110 fn derivatives(&self, state: [f64; 4], current: f64) -> [f64; 4] {
2111 let [voltage, h_gate, n_gate, s_gate] = state;
2112 if !Self::finite_values(&state)
2113 || !Self::gate_in_range(h_gate)
2114 || !Self::gate_in_range(n_gate)
2115 || !Self::gate_in_range(s_gate)
2116 {
2117 return [f64::NAN; 4];
2118 }
2119 let rates = self.rates(voltage);
2120 let i_na = self.g_na * rates[0].powi(3) * h_gate * (voltage - self.e_na);
2121 let i_k = self.g_k * n_gate.powi(4) * (voltage - self.e_k);
2122 let i_s = self.g_s * s_gate * (voltage - self.e_s);
2123 let i_l = self.g_l * (voltage - self.e_l);
2124 [
2125 -i_na - i_k - i_s - i_l + current,
2126 (rates[1] - h_gate) / rates[4],
2127 (rates[2] - n_gate) / rates[5],
2128 (rates[3] - s_gate) / rates[6],
2129 ]
2130 }
2131
2132 fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
2133 [
2134 state[0] + scale * slope[0],
2135 state[1] + scale * slope[1],
2136 state[2] + scale * slope[2],
2137 state[3] + scale * slope[3],
2138 ]
2139 }
2140
2141 fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
2142 let state = [self.v, self.h, self.n, self.s];
2143 let half_dt = 0.5 * self.dt;
2144 let k1 = self.derivatives(state, current);
2145 let k2 = self.derivatives(Self::add_scaled(state, k1, half_dt), current);
2146 let k3 = self.derivatives(Self::add_scaled(state, k2, half_dt), current);
2147 let k4 = self.derivatives(Self::add_scaled(state, k3, self.dt), current);
2148 let candidate = [
2149 state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
2150 state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
2151 state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
2152 state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
2153 ];
2154 if Self::finite_values(&candidate)
2155 && Self::gate_in_range(candidate[1])
2156 && Self::gate_in_range(candidate[2])
2157 && Self::gate_in_range(candidate[3])
2158 {
2159 Some(candidate)
2160 } else {
2161 None
2162 }
2163 }
2164
2165 pub fn step(&mut self, current: f64) -> i32 {
2166 if !current.is_finite() || !self.valid_runtime() {
2167 return 0;
2168 }
2169 let v_prev = self.v;
2170 let Some(candidate) = self.rk4_candidate(current) else {
2171 return 0;
2172 };
2173 self.v = candidate[0];
2174 self.h = candidate[1];
2175 self.n = candidate[2];
2176 self.s = candidate[3];
2177 if self.v >= self.v_threshold && v_prev < self.v_threshold {
2178 1
2179 } else {
2180 0
2181 }
2182 }
2183
2184 pub fn reset(&mut self) {
2185 self.v = -60.0;
2186 self.h = 0.6;
2187 self.n = 0.3;
2188 self.s = 0.5;
2189 }
2190}
2191impl Default for AvRonCardiacNeuron {
2192 fn default() -> Self {
2193 Self::new()
2194 }
2195}
2196
2197#[derive(Clone, Debug)]
2199pub struct DurstewitzDopamineNeuron {
2200 pub v: f64,
2201 pub h_na: f64,
2202 pub n_k: f64,
2203 pub g_na: f64,
2204 pub g_k: f64,
2205 pub g_nmda: f64,
2206 pub g_l: f64,
2207 pub e_na: f64,
2208 pub e_k: f64,
2209 pub e_nmda: f64,
2210 pub e_l: f64,
2211 pub mg: f64,
2212 pub d1_level: f64,
2213 pub g_nmda_scale: f64,
2214 pub g_k_scale: f64,
2215 pub v_shift_na: f64,
2216 pub dt: f64,
2217 pub v_threshold: f64,
2218}
2219
2220impl DurstewitzDopamineNeuron {
2221 pub fn new() -> Self {
2222 Self {
2223 v: -65.0,
2224 h_na: 0.7,
2225 n_k: 0.2,
2226 g_na: 45.0,
2227 g_k: 18.0,
2228 g_nmda: 0.5,
2229 g_l: 0.02,
2230 e_na: 55.0,
2231 e_k: -80.0,
2232 e_nmda: 0.0,
2233 e_l: -65.0,
2234 mg: 1.0,
2235 d1_level: 0.0,
2236 g_nmda_scale: 2.5,
2237 g_k_scale: 1.5,
2238 v_shift_na: -5.0,
2239 dt: 0.05,
2240 v_threshold: -20.0,
2241 }
2242 }
2243 fn derivatives(&self, v: f64, h_na: f64, n_k: f64, current: f64) -> [f64; 3] {
2250 let v_sh = self.d1_level * self.v_shift_na;
2251 let m_na_inf = 1.0 / (1.0 + (-(v + 30.0 + v_sh) / 9.5).exp());
2252 let h_na_inf = 1.0 / (1.0 + ((v + 53.0) / 7.0).exp());
2253 let n_k_inf = 1.0 / (1.0 + (-(v + 30.0) / 10.0).exp());
2254 let tau_h = 0.5 + 14.0 / (1.0 + ((v + 50.0) / 12.0).exp());
2255 let tau_n = 1.0 + 11.0 / (1.0 + ((v + 40.0) / 10.0).exp());
2256 let d_h_na = (h_na_inf - h_na) / tau_h;
2257 let d_n_k = (n_k_inf - n_k) / tau_n;
2258 let mg_block = 1.0 / (1.0 + self.mg / 3.57 * (-0.062 * v).exp());
2259 let nmda_g = self.g_nmda * (1.0 + self.d1_level * (self.g_nmda_scale - 1.0));
2260 let k_g = self.g_k * (1.0 + self.d1_level * (self.g_k_scale - 1.0));
2261 let i_na = self.g_na * m_na_inf * m_na_inf * m_na_inf * h_na * (v - self.e_na);
2262 let i_k = k_g * n_k * n_k * n_k * n_k * (v - self.e_k);
2263 let i_nmda = nmda_g * mg_block * (v - self.e_nmda);
2264 let i_l = self.g_l * (v - self.e_l);
2265 let d_v = -i_na - i_k - i_nmda - i_l + current;
2266 [d_v, d_h_na, d_n_k]
2267 }
2268
2269 fn rk4_substep(&self, s: [f64; 3], current: f64) -> [f64; 3] {
2271 let dt = self.dt;
2272 let k1 = self.derivatives(s[0], s[1], s[2], current);
2273 let k2 = self.derivatives(
2274 s[0] + 0.5 * dt * k1[0],
2275 s[1] + 0.5 * dt * k1[1],
2276 s[2] + 0.5 * dt * k1[2],
2277 current,
2278 );
2279 let k3 = self.derivatives(
2280 s[0] + 0.5 * dt * k2[0],
2281 s[1] + 0.5 * dt * k2[1],
2282 s[2] + 0.5 * dt * k2[2],
2283 current,
2284 );
2285 let k4 = self.derivatives(
2286 s[0] + dt * k3[0],
2287 s[1] + dt * k3[1],
2288 s[2] + dt * k3[2],
2289 current,
2290 );
2291 let mut out = [0.0_f64; 3];
2292 for i in 0..3 {
2293 out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
2294 }
2295 out
2296 }
2297
2298 pub fn step(&mut self, current: f64) -> i32 {
2299 let v_prev = self.v;
2300 let s = self.rk4_substep([self.v, self.h_na, self.n_k], current);
2301 self.v = s[0];
2302 self.h_na = s[1];
2303 self.n_k = s[2];
2304 if self.v >= self.v_threshold && v_prev < self.v_threshold {
2305 1
2306 } else {
2307 0
2308 }
2309 }
2310 pub fn reset(&mut self) {
2311 self.v = -65.0;
2312 self.h_na = 0.7;
2313 self.n_k = 0.2;
2314 }
2315}
2316impl Default for DurstewitzDopamineNeuron {
2317 fn default() -> Self {
2318 Self::new()
2319 }
2320}
2321
2322#[derive(Clone, Debug)]
2324pub struct HillTononiNeuron {
2325 pub v: f64,
2326 pub h_na: f64,
2327 pub n_k: f64,
2328 pub m_h: f64,
2329 pub h_t: f64,
2330 pub na_i: f64,
2331 pub g_na: f64,
2332 pub g_k: f64,
2333 pub g_h: f64,
2334 pub g_t: f64,
2335 pub g_kna: f64,
2336 pub g_l: f64,
2337 pub e_na: f64,
2338 pub e_k: f64,
2339 pub e_h: f64,
2340 pub e_ca: f64,
2341 pub e_l: f64,
2342 pub na_pump_max: f64,
2343 pub na_eq: f64,
2344 pub dt: f64,
2345 pub v_threshold: f64,
2346}
2347
2348impl HillTononiNeuron {
2349 pub fn new() -> Self {
2350 Self {
2351 v: -65.0,
2352 h_na: 0.6,
2353 n_k: 0.3,
2354 m_h: 0.0,
2355 h_t: 0.9,
2356 na_i: 5.0,
2357 g_na: 50.0,
2358 g_k: 5.0,
2359 g_h: 1.0,
2360 g_t: 3.0,
2361 g_kna: 1.33,
2362 g_l: 0.02,
2363 e_na: 50.0,
2364 e_k: -90.0,
2365 e_h: -43.0,
2366 e_ca: 120.0,
2367 e_l: -70.0,
2368 na_pump_max: 20.0,
2369 na_eq: 9.5,
2370 dt: 0.05,
2371 v_threshold: -20.0,
2372 }
2373 }
2374 fn derivatives(
2381 &self,
2382 v: f64,
2383 h_na: f64,
2384 n_k: f64,
2385 m_h: f64,
2386 h_t: f64,
2387 na_i: f64,
2388 current: f64,
2389 ) -> [f64; 6] {
2390 let m_na_inf = 1.0 / (1.0 + (-(v + 38.0) / 10.0).exp());
2391 let h_na_inf = 1.0 / (1.0 + ((v + 43.0) / 6.0).exp());
2392 let n_k_inf = 1.0 / (1.0 + (-(v + 27.0) / 11.5).exp());
2393 let m_h_inf = 1.0 / (1.0 + ((v + 75.0) / 5.5).exp());
2394 let m_t_inf = 1.0 / (1.0 + (-(v + 59.0) / 6.2).exp());
2395 let h_t_inf = 1.0 / (1.0 + ((v + 83.0) / 4.0).exp());
2396 let hill_base = 38.7 / na_i.max(0.01);
2397 let w_kna = 0.37 / (1.0 + hill_base * hill_base * hill_base * hill_base.sqrt());
2398 let tau_h_na = (1.0 + 10.0 / (1.0 + ((v + 40.0) / 10.0).exp())).max(0.1);
2399 let z_n = (v + 50.0) / 25.0;
2400 let tau_n_k = (5.0 + 47.0 * (-(z_n * z_n)).exp()).max(0.1);
2401 let tau_m_h =
2402 (20.0 + 1000.0 / (((v + 71.5) / 14.2).exp() + (-(v + 89.0) / 11.6).exp())).max(1.0);
2403 let tau_h_t = if v < -81.0 {
2404 (30.8 + 211.4 * ((v + 115.2) / 5.0).exp() / (1.0 + ((v + 86.0) / 3.2).exp())).max(0.1)
2405 } else {
2406 10.0
2407 };
2408 let d_h_na = (h_na_inf - h_na) / tau_h_na;
2409 let d_n_k = (n_k_inf - n_k) / tau_n_k;
2410 let d_m_h = (m_h_inf - m_h) / tau_m_h;
2411 let d_h_t = (h_t_inf - h_t) / tau_h_t;
2412 let i_na = self.g_na * m_na_inf * m_na_inf * m_na_inf * h_na * (v - self.e_na);
2413 let i_k = self.g_k * n_k * n_k * n_k * n_k * (v - self.e_k);
2414 let i_h = self.g_h * m_h * (v - self.e_h);
2415 let i_t = self.g_t * m_t_inf * m_t_inf * h_t * (v - self.e_ca);
2416 let i_kna = self.g_kna * w_kna * (v - self.e_k);
2417 let i_l = self.g_l * (v - self.e_l);
2418 let d_v = -i_na - i_k - i_h - i_t - i_kna - i_l + current;
2419 let d_na_i = -0.001 * i_na - self.na_pump_max * (na_i / (na_i + self.na_eq));
2420 [d_v, d_h_na, d_n_k, d_m_h, d_h_t, d_na_i]
2421 }
2422
2423 fn rk4_substep(&self, s: [f64; 6], current: f64) -> [f64; 6] {
2425 let dt = self.dt;
2426 let k1 = self.derivatives(s[0], s[1], s[2], s[3], s[4], s[5], current);
2427 let k2 = self.derivatives(
2428 s[0] + 0.5 * dt * k1[0],
2429 s[1] + 0.5 * dt * k1[1],
2430 s[2] + 0.5 * dt * k1[2],
2431 s[3] + 0.5 * dt * k1[3],
2432 s[4] + 0.5 * dt * k1[4],
2433 s[5] + 0.5 * dt * k1[5],
2434 current,
2435 );
2436 let k3 = self.derivatives(
2437 s[0] + 0.5 * dt * k2[0],
2438 s[1] + 0.5 * dt * k2[1],
2439 s[2] + 0.5 * dt * k2[2],
2440 s[3] + 0.5 * dt * k2[3],
2441 s[4] + 0.5 * dt * k2[4],
2442 s[5] + 0.5 * dt * k2[5],
2443 current,
2444 );
2445 let k4 = self.derivatives(
2446 s[0] + dt * k3[0],
2447 s[1] + dt * k3[1],
2448 s[2] + dt * k3[2],
2449 s[3] + dt * k3[3],
2450 s[4] + dt * k3[4],
2451 s[5] + dt * k3[5],
2452 current,
2453 );
2454 let mut out = [0.0_f64; 6];
2455 for i in 0..6 {
2456 out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
2457 }
2458 out
2459 }
2460
2461 pub fn step(&mut self, current: f64) -> i32 {
2462 let v_prev = self.v;
2463 let s = self.rk4_substep(
2464 [self.v, self.h_na, self.n_k, self.m_h, self.h_t, self.na_i],
2465 current,
2466 );
2467 self.v = s[0];
2468 self.h_na = s[1];
2469 self.n_k = s[2];
2470 self.m_h = s[3];
2471 self.h_t = s[4];
2472 self.na_i = s[5].max(0.0);
2473 if self.v >= self.v_threshold && v_prev < self.v_threshold {
2474 1
2475 } else {
2476 0
2477 }
2478 }
2479 pub fn reset(&mut self) {
2480 self.v = -65.0;
2481 self.h_na = 0.6;
2482 self.n_k = 0.3;
2483 self.m_h = 0.0;
2484 self.h_t = 0.9;
2485 self.na_i = 5.0;
2486 }
2487}
2488impl Default for HillTononiNeuron {
2489 fn default() -> Self {
2490 Self::new()
2491 }
2492}
2493
2494#[derive(Clone, Debug)]
2496pub struct BertramPhantomBurster {
2497 pub v: f64,
2498 pub s1: f64,
2499 pub s2: f64,
2500 pub g_ca: f64,
2501 pub g_k: f64,
2502 pub g_s1: f64,
2503 pub g_s2: f64,
2504 pub g_l: f64,
2505 pub e_ca: f64,
2506 pub e_k: f64,
2507 pub e_l: f64,
2508 pub c_m: f64,
2509 pub v_m: f64,
2510 pub s_m: f64,
2511 pub v_n: f64,
2512 pub s_n: f64,
2513 pub v_s1: f64,
2514 pub s_s1: f64,
2515 pub v_s2: f64,
2516 pub s_s2: f64,
2517 pub tau_s1: f64,
2518 pub tau_s2: f64,
2519 pub dt: f64,
2520 pub v_threshold: f64,
2521}
2522
2523impl BertramPhantomBurster {
2524 pub fn new() -> Self {
2525 Self {
2526 v: -50.0,
2527 s1: 0.1,
2528 s2: 0.1,
2529 g_ca: 3.6,
2530 g_k: 10.0,
2531 g_s1: 4.0,
2532 g_s2: 4.0,
2533 g_l: 0.2,
2534 e_ca: 25.0,
2535 e_k: -75.0,
2536 e_l: -40.0,
2537 c_m: 5.3,
2538 v_m: -20.0,
2539 s_m: 12.0,
2540 v_n: -16.0,
2541 s_n: 5.6,
2542 v_s1: -40.0,
2543 s_s1: 10.0,
2544 v_s2: -42.0,
2545 s_s2: 0.4,
2546 tau_s1: 20000.0,
2547 tau_s2: 100000.0,
2548 dt: 0.5,
2549 v_threshold: -20.0,
2550 }
2551 }
2552 pub fn step(&mut self, current: f64) -> i32 {
2553 let v_prev = self.v;
2554 let m_inf = 1.0 / (1.0 + (-(self.v - self.v_m) / self.s_m).exp());
2555 let n_inf = 1.0 / (1.0 + (-(self.v - self.v_n) / self.s_n).exp());
2556 let s1_inf = 1.0 / (1.0 + (-(self.v - self.v_s1) / self.s_s1).exp());
2557 let s2_inf = 1.0 / (1.0 + (-(self.v - self.v_s2) / self.s_s2).exp());
2558 let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
2559 let i_k = self.g_k * n_inf * (self.v - self.e_k);
2560 let i_s1 = self.g_s1 * self.s1 * (self.v - self.e_k);
2561 let i_s2 = self.g_s2 * self.s2 * (self.v - self.e_k);
2562 let i_l = self.g_l * (self.v - self.e_l);
2563 self.v += (-i_ca - i_k - i_s1 - i_s2 - i_l + current) / self.c_m * self.dt;
2564 self.s1 += (s1_inf - self.s1) / self.tau_s1 * self.dt;
2565 self.s2 += (s2_inf - self.s2) / self.tau_s2 * self.dt;
2566 if self.v >= self.v_threshold && v_prev < self.v_threshold {
2567 1
2568 } else {
2569 0
2570 }
2571 }
2572 pub fn reset(&mut self) {
2573 self.v = -50.0;
2574 self.s1 = 0.1;
2575 self.s2 = 0.1;
2576 }
2577}
2578impl Default for BertramPhantomBurster {
2579 fn default() -> Self {
2580 Self::new()
2581 }
2582}
2583
2584#[derive(Clone, Debug)]
2586pub struct YamadaNeuron {
2587 pub v: f64,
2588 pub n: f64,
2589 pub q: f64,
2590 pub g_na: f64,
2591 pub g_k: f64,
2592 pub g_q: f64,
2593 pub g_l: f64,
2594 pub e_na: f64,
2595 pub e_k: f64,
2596 pub e_q: f64,
2597 pub e_l: f64,
2598 pub tau_q: f64,
2599 pub dt: f64,
2600 pub v_threshold: f64,
2601}
2602
2603impl YamadaNeuron {
2604 pub fn new() -> Self {
2605 Self {
2606 v: -60.0,
2607 n: 0.1,
2608 q: 0.0,
2609 g_na: 20.0,
2610 g_k: 10.0,
2611 g_q: 5.0,
2612 g_l: 0.5,
2613 e_na: 60.0,
2614 e_k: -80.0,
2615 e_q: -80.0,
2616 e_l: -60.0,
2617 tau_q: 300.0,
2618 dt: 0.05,
2619 v_threshold: -20.0,
2620 }
2621 }
2622 pub fn step(&mut self, current: f64) -> i32 {
2623 let v_prev = self.v;
2624 let m_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 9.5).exp());
2625 let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 10.0).exp());
2626 let q_inf = 1.0 / (1.0 + (-(self.v + 50.0) / 10.0).exp());
2627 let tau_n = 1.0 + 7.5 / (1.0 + ((self.v + 40.0) / 12.0).exp());
2628 let i_na = self.g_na * m_inf.powi(3) * (1.0 - self.n) * (self.v - self.e_na);
2629 let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
2630 let i_q = self.g_q * self.q * (self.v - self.e_q);
2631 let i_l = self.g_l * (self.v - self.e_l);
2632 self.v += (-i_na - i_k - i_q - i_l + current) * self.dt;
2633 self.n += (n_inf - self.n) / tau_n * self.dt;
2634 self.q += (q_inf - self.q) / self.tau_q * self.dt;
2635 if self.v >= self.v_threshold && v_prev < self.v_threshold {
2636 1
2637 } else {
2638 0
2639 }
2640 }
2641 pub fn reset(&mut self) {
2642 self.v = -60.0;
2643 self.n = 0.1;
2644 self.q = 0.0;
2645 }
2646}
2647impl Default for YamadaNeuron {
2648 fn default() -> Self {
2649 Self::new()
2650 }
2651}
2652
2653#[cfg(test)]
2654mod tests {
2655 use super::*;
2656
2657 #[test]
2658 fn hh_fires() {
2659 let mut n = HodgkinHuxleyNeuron::new();
2660 let t: i32 = (0..100).map(|_| n.step(10.0)).sum();
2661 assert!(t > 0);
2662 }
2663 #[test]
2664 fn traub_fires() {
2665 let mut n = TraubMilesNeuron::new();
2666 let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
2667 assert!(t > 0);
2668 }
2669 #[test]
2670 fn wb_fires() {
2671 let mut n = WangBuzsakiNeuron::new();
2672 let t: i32 = (0..200).map(|_| n.step(2.0)).sum();
2673 assert!(t > 0);
2674 }
2675 #[test]
2676 fn cs_fires() {
2677 let mut n = ConnorStevensNeuron::new();
2678 let t: i32 = (0..200).map(|_| n.step(10.0)).sum();
2679 assert!(t > 0);
2680 }
2681 #[test]
2682 fn destexhe_fires() {
2683 let mut n = DestexheThalamicNeuron::new();
2684 let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
2685 assert!(t > 0);
2686 }
2687 #[test]
2688 fn hb_fires() {
2689 let mut n = HuberBraunNeuron::new();
2690 let t: i32 = (0..5000).map(|_| n.step(10.0)).sum();
2691 assert!(t > 0);
2692 }
2693 #[test]
2694 fn golomb_fires() {
2695 let mut n = GolombFSNeuron::new();
2696 let t: i32 = (0..2000).map(|_| n.step(200.0)).sum();
2697 assert!(t > 0);
2698 }
2699 #[test]
2700 fn pospischil_fires() {
2701 let mut n = PospischilNeuron::new();
2702 let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
2703 assert!(t > 0);
2704 }
2705 #[test]
2706 fn mainen_fires() {
2707 let mut n = MainenSejnowskiNeuron::new();
2708 let t: i32 = (0..5000).map(|_| n.step(500.0)).sum();
2709 assert!(t > 0);
2710 }
2711 #[test]
2712 fn purkinje_fires() {
2713 let mut n = DeSchutterPurkinjeNeuron::new();
2714 let t: i32 = (0..5000).map(|_| n.step(200.0)).sum();
2715 assert!(t > 0);
2716 }
2717 #[test]
2718 fn plant_r15_fires() {
2719 let mut n = PlantR15Neuron::new();
2720 let t: i32 = (0..500).map(|_| n.step(2.0)).sum();
2721 assert!(t > 0);
2722 }
2723 #[test]
2724 fn prescott_fires() {
2725 let mut n = PrescottNeuron::new();
2726 let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
2727 assert!(t > 0);
2728 }
2729 #[test]
2730 fn mn_fires() {
2731 let mut n = MihalasNieburNeuron::new();
2732 let t: i32 = (0..100).map(|_| n.step(5.0)).sum();
2733 assert!(t > 0);
2734 }
2735 #[test]
2736 fn glif_fires() {
2737 let mut n = GLIFNeuron::new();
2738 let t: i32 = (0..200).map(|_| n.step(30.0)).sum();
2739 assert!(t > 0);
2740 }
2741 #[test]
2742 fn gif_pop_fires() {
2743 let mut n = GIFPopulationNeuron::new(42);
2744 let t: i32 = (0..1000).map(|_| n.step(30.0)).sum();
2745 assert!(t > 0);
2746 }
2747 #[test]
2748 fn avron_fires() {
2749 let mut n = AvRonCardiacNeuron::new();
2750 let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
2751 assert!(t > 0);
2752 }
2753 #[test]
2754 fn durstewitz_fires() {
2755 let mut n = DurstewitzDopamineNeuron::new();
2756 let t: i32 = (0..1000).map(|_| n.step(3.0)).sum();
2757 assert!(t > 0);
2758 }
2759 #[test]
2760 fn hill_tononi_fires() {
2761 let mut n = HillTononiNeuron::new();
2762 let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
2763 assert!(t > 0);
2764 }
2765 #[test]
2766 fn bertram_fires() {
2767 let mut n = BertramPhantomBurster::new();
2768 let t: i32 = (0..10000).map(|_| n.step(200.0)).sum();
2769 assert!(t > 0);
2770 }
2771 #[test]
2772 fn yamada_fires() {
2773 let mut n = YamadaNeuron::new();
2774 let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
2775 assert!(t > 0);
2776 }
2777
2778 #[test]
2782 fn hh_silent_without_input() {
2783 let mut n = HodgkinHuxleyNeuron::new();
2784 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
2785 assert_eq!(t, 0);
2786 }
2787 #[test]
2788 fn hh_reset_clears_state() {
2789 let mut n = HodgkinHuxleyNeuron::new();
2790 for _ in 0..100 {
2791 n.step(10.0);
2792 }
2793 n.reset();
2794 assert!((n.v - (-65.0)).abs() < 1e-10);
2795 assert!((n.m - 0.05).abs() < 1e-10);
2796 assert!((n.h - 0.6).abs() < 1e-10);
2797 assert!((n.n - 0.32).abs() < 1e-10);
2798 }
2799 #[test]
2800 fn hh_extreme_input_bounded() {
2801 let mut n = HodgkinHuxleyNeuron::new();
2802 for _ in 0..200 {
2803 n.step(1e4);
2804 }
2805 assert!(n.v.is_finite());
2806 }
2807 #[test]
2808 fn hh_gates_bounded() {
2809 let mut n = HodgkinHuxleyNeuron::new();
2810 for _ in 0..500 {
2811 n.step(10.0);
2812 }
2813 assert!(n.m >= 0.0 && n.m <= 1.0, "m={}", n.m);
2814 assert!(n.h >= 0.0 && n.h <= 1.0, "h={}", n.h);
2815 assert!(n.n >= 0.0 && n.n <= 1.0, "n={}", n.n);
2816 }
2817 #[test]
2818 fn hh_negative_input_no_crash() {
2819 let mut n = HodgkinHuxleyNeuron::new();
2820 for _ in 0..200 {
2821 n.step(-20.0);
2822 }
2823 assert!(n.v.is_finite());
2824 }
2825 #[test]
2826 fn hh_nan_input_no_panic() {
2827 let mut n = HodgkinHuxleyNeuron::new();
2828 n.step(f64::NAN);
2829 }
2830 #[test]
2831 fn hh_sodium_potassium_opposition() {
2832 let mut n = HodgkinHuxleyNeuron::new();
2834 for _ in 0..50 {
2835 n.step(10.0);
2836 }
2837 assert!(n.n > 0.32, "K activation n should increase during spiking");
2839 }
2840
2841 #[test]
2843 fn traub_silent_without_input() {
2844 let mut n = TraubMilesNeuron::new();
2845 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
2846 assert_eq!(t, 0);
2847 }
2848 #[test]
2849 fn traub_reset_clears_state() {
2850 let mut n = TraubMilesNeuron::new();
2851 for _ in 0..100 {
2852 n.step(5.0);
2853 }
2854 n.reset();
2855 assert!((n.v - (-67.0)).abs() < 1e-10);
2856 }
2857 #[test]
2858 fn traub_extreme_bounded() {
2859 let mut n = TraubMilesNeuron::new();
2860 for _ in 0..200 {
2861 n.step(1e4);
2862 }
2863 assert!(n.v.is_finite());
2864 }
2865 #[test]
2866 fn traub_gates_bounded() {
2867 let mut n = TraubMilesNeuron::new();
2868 for _ in 0..500 {
2869 n.step(5.0);
2870 }
2871 assert!(n.m >= 0.0 && n.m <= 1.01);
2872 assert!(n.h >= 0.0 && n.h <= 1.01);
2873 assert!(n.n >= 0.0 && n.n <= 1.01);
2874 }
2875 #[test]
2876 fn traub_weak_negative_no_crash() {
2877 let mut n = TraubMilesNeuron::new();
2878 for _ in 0..200 {
2879 n.step(-5.0);
2880 }
2881 assert!(n.v.is_finite());
2882 }
2883 #[test]
2884 fn traub_nan_no_panic() {
2885 let mut n = TraubMilesNeuron::new();
2886 n.step(f64::NAN);
2887 }
2888 #[test]
2889 fn traub_rk4_reference_point() {
2890 let mut n = TraubMilesNeuron::new();
2891 n.v = -63.5;
2892 n.m = 0.08;
2893 n.h = 0.55;
2894 n.n = 0.32;
2895 let spike = n.step(4.0);
2896 assert_eq!(spike, 0);
2897 assert!((n.v - (-65.6638958700765)).abs() < 1e-13);
2898 assert!((n.m - 0.04237301812907925).abs() < 1e-15);
2899 assert!((n.h - 0.5626824931070477).abs() < 1e-15);
2900 assert!((n.n - 0.30356298261126924).abs() < 1e-15);
2901 assert!((n.v - (-65.66233161606698)).abs() > 1e-3);
2902 }
2903
2904 #[test]
2906 fn wb_silent_without_input() {
2907 let mut n = WangBuzsakiNeuron::new();
2908 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
2909 assert_eq!(t, 0);
2910 }
2911 #[test]
2912 fn wb_reset_clears_state() {
2913 let mut n = WangBuzsakiNeuron::new();
2914 for _ in 0..100 {
2915 n.step(2.0);
2916 }
2917 n.reset();
2918 assert!((n.v - (-65.0)).abs() < 1e-10);
2919 }
2920 #[test]
2921 fn wb_extreme_bounded() {
2922 let mut n = WangBuzsakiNeuron::new();
2923 for _ in 0..200 {
2924 n.step(1e4);
2925 }
2926 assert!(n.v.is_finite());
2927 }
2928 #[test]
2929 fn wb_fast_spiking_high_rate() {
2930 let mut n = WangBuzsakiNeuron::new();
2932 let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
2933 assert!(t > 10, "WB FS should produce many spikes, got {}", t);
2934 }
2935 #[test]
2936 fn wb_gates_bounded() {
2937 let mut n = WangBuzsakiNeuron::new();
2938 for _ in 0..500 {
2939 n.step(2.0);
2940 }
2941 assert!(n.h >= 0.0 && n.h <= 1.0);
2942 assert!(n.n >= 0.0 && n.n <= 1.0);
2943 }
2944 #[test]
2945 fn wb_negative_no_crash() {
2946 let mut n = WangBuzsakiNeuron::new();
2947 for _ in 0..200 {
2948 n.step(-10.0);
2949 }
2950 assert!(n.v.is_finite());
2951 }
2952 #[test]
2953 fn wb_nan_no_panic() {
2954 let mut n = WangBuzsakiNeuron::new();
2955 n.step(f64::NAN);
2956 }
2957
2958 #[test]
2960 fn cs_silent_without_input() {
2961 let mut n = ConnorStevensNeuron::new();
2962 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
2963 assert_eq!(t, 0);
2964 }
2965 #[test]
2966 fn cs_reset_clears_state() {
2967 let mut n = ConnorStevensNeuron::new();
2968 for _ in 0..100 {
2969 n.step(10.0);
2970 }
2971 n.reset();
2972 assert!((n.v - (-68.0)).abs() < 1e-10);
2973 }
2974 #[test]
2975 fn cs_extreme_bounded() {
2976 let mut n = ConnorStevensNeuron::new();
2977 for _ in 0..200 {
2978 n.step(1e4);
2979 }
2980 assert!(n.v.is_finite());
2981 }
2982 #[test]
2983 fn cs_a_type_delays_spike() {
2984 let mut n_with_a = ConnorStevensNeuron::new();
2987 let mut n_no_a = ConnorStevensNeuron::new();
2988 n_no_a.g_a = 0.0;
2989 let spikes_with: i32 = (0..50).map(|_| n_with_a.step(8.0)).sum();
2990 let spikes_no: i32 = (0..50).map(|_| n_no_a.step(8.0)).sum();
2991 assert!(
2993 spikes_no >= spikes_with,
2994 "without A-type K, should fire more: no_a={} vs with_a={}",
2995 spikes_no,
2996 spikes_with
2997 );
2998 }
2999 #[test]
3000 fn cs_gates_bounded() {
3001 let mut n = ConnorStevensNeuron::new();
3002 for _ in 0..500 {
3003 n.step(10.0);
3004 }
3005 assert!(n.a >= 0.0 && n.a <= 1.5, "a={}", n.a); assert!(n.b >= 0.0 && n.b <= 1.0, "b={}", n.b);
3007 }
3008 #[test]
3009 fn cs_negative_no_crash() {
3010 let mut n = ConnorStevensNeuron::new();
3011 for _ in 0..200 {
3012 n.step(-20.0);
3013 }
3014 assert!(n.v.is_finite());
3015 }
3016 #[test]
3017 fn cs_invalid_current_preserves_state() {
3018 let mut n = ConnorStevensNeuron::new();
3019 let before = (n.v, n.m, n.h, n.n, n.a, n.b);
3020 assert_eq!(n.step(f64::NAN), 0);
3021 assert_eq!((n.v, n.m, n.h, n.n, n.a, n.b), before);
3022 }
3023 #[test]
3024 fn cs_corrupt_runtime_state_preserves_state() {
3025 let mut n = ConnorStevensNeuron::new();
3026 n.b = f64::INFINITY;
3027 let before = (n.v, n.m, n.h, n.n, n.a, n.b);
3028 assert_eq!(n.step(6.0), 0);
3029 assert_eq!((n.v, n.m, n.h, n.n, n.a, n.b), before);
3030 }
3031
3032 #[test]
3034 fn destexhe_no_crash_zero_input() {
3035 let mut n = DestexheThalamicNeuron::new();
3036 let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3037 assert!(n.v.is_finite());
3039 }
3040 #[test]
3041 fn destexhe_reset_clears_state() {
3042 let mut n = DestexheThalamicNeuron::new();
3043 for _ in 0..200 {
3044 n.step(5.0);
3045 }
3046 n.reset();
3047 assert!((n.v - (-65.0)).abs() < 1e-10);
3048 }
3049 #[test]
3050 fn destexhe_extreme_bounded() {
3051 let mut n = DestexheThalamicNeuron::new();
3052 for _ in 0..200 {
3053 n.step(1e4);
3054 }
3055 assert!(n.v.is_finite());
3056 }
3057 #[test]
3058 fn destexhe_t_current_rebound() {
3059 let v_hyp = -90.0_f64;
3063 let h_t_inf = 1.0 / (1.0 + ((v_hyp + 81.0) / 4.0).exp());
3064 assert!(h_t_inf > 0.9, "h_t_inf should be ~1 at V=-90: {}", h_t_inf);
3065
3066 let mut n = DestexheThalamicNeuron::new();
3068 for _ in 0..500 {
3069 n.step(-5.0);
3070 }
3071 let v_before_release = n.v;
3072 for _ in 0..500 {
3073 n.step(0.0);
3074 }
3075 assert!(n.v.is_finite());
3076 assert!(
3077 n.v > v_before_release,
3078 "V should increase after release (rebound)"
3079 );
3080 }
3081 #[test]
3082 fn destexhe_negative_no_crash() {
3083 let mut n = DestexheThalamicNeuron::new();
3084 for _ in 0..200 {
3085 n.step(-20.0);
3086 }
3087 assert!(n.v.is_finite());
3088 }
3089 #[test]
3090 fn destexhe_nan_no_panic() {
3091 let mut n = DestexheThalamicNeuron::new();
3092 n.step(f64::NAN);
3093 }
3094
3095 #[test]
3097 fn hb_silent_without_input() {
3098 let mut n = HuberBraunNeuron::new();
3099 let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3100 assert!(n.v.is_finite());
3102 }
3103 #[test]
3104 fn hb_reset_clears_state() {
3105 let mut n = HuberBraunNeuron::new();
3106 for _ in 0..200 {
3107 n.step(10.0);
3108 }
3109 n.reset();
3110 assert!((n.v - (-50.0)).abs() < 1e-10);
3111 }
3112 #[test]
3113 fn hb_extreme_bounded() {
3114 let mut n = HuberBraunNeuron::new();
3115 for _ in 0..200 {
3116 n.step(1e4);
3117 }
3118 assert!(n.v.is_finite());
3119 }
3120 #[test]
3121 fn hb_negative_no_crash() {
3122 let mut n = HuberBraunNeuron::new();
3123 for _ in 0..500 {
3124 n.step(-10.0);
3125 }
3126 assert!(n.v.is_finite());
3127 }
3128 #[test]
3129 fn hb_nan_no_panic() {
3130 let mut n = HuberBraunNeuron::new();
3131 n.step(f64::NAN);
3132 }
3133
3134 #[test]
3136 fn golomb_silent_without_input() {
3137 let mut n = GolombFSNeuron::new();
3138 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
3139 assert_eq!(t, 0);
3140 }
3141 #[test]
3142 fn golomb_reset_clears_state() {
3143 let mut n = GolombFSNeuron::new();
3144 for _ in 0..100 {
3145 n.step(200.0);
3146 }
3147 n.reset();
3148 assert!((n.v - (-65.0)).abs() < 1e-10);
3149 }
3150 #[test]
3151 fn golomb_extreme_bounded() {
3152 let mut n = GolombFSNeuron::new();
3154 for _ in 0..200 {
3155 n.step(200.0);
3156 }
3157 assert!(n.v.is_finite());
3158 }
3159 #[test]
3160 fn golomb_kv3_enables_fast_spiking() {
3161 let mut n = GolombFSNeuron::new();
3163 let t: i32 = (0..5000).map(|_| n.step(300.0)).sum();
3164 assert!(t > 0, "Golomb FS should fire with high input, got {}", t);
3165 }
3166 #[test]
3167 fn golomb_negative_no_crash() {
3168 let mut n = GolombFSNeuron::new();
3169 for _ in 0..200 {
3170 n.step(-100.0);
3171 }
3172 assert!(n.v.is_finite());
3173 }
3174 #[test]
3175 fn golomb_nan_no_panic() {
3176 let mut n = GolombFSNeuron::new();
3177 n.step(f64::NAN);
3178 }
3179
3180 #[test]
3182 fn pospischil_silent_without_input() {
3183 let mut n = PospischilNeuron::new();
3184 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
3185 assert_eq!(t, 0);
3186 }
3187 #[test]
3188 fn pospischil_reset_clears_state() {
3189 let mut n = PospischilNeuron::new();
3190 for _ in 0..100 {
3191 n.step(5.0);
3192 }
3193 n.reset();
3194 assert!((n.v - (-70.0)).abs() < 1e-10);
3195 }
3196 #[test]
3197 fn pospischil_moderate_input_stable() {
3198 let mut n = PospischilNeuron::new();
3199 for _ in 0..200 {
3200 n.step(10.0);
3201 }
3202 assert!(n.v.is_finite());
3203 }
3204 #[test]
3205 fn pospischil_m_current_present() {
3206 let mut n = PospischilNeuron::new();
3207 for _ in 0..200 {
3208 n.step(5.0);
3209 }
3210 assert!(n.p > 0.0, "M-current (p) should activate during spiking");
3211 }
3212 #[test]
3213 fn pospischil_negative_no_crash() {
3214 let mut n = PospischilNeuron::new();
3215 for _ in 0..200 {
3216 n.step(-10.0);
3217 }
3218 assert!(n.v.is_finite());
3219 }
3220 #[test]
3221 fn pospischil_nan_no_panic() {
3222 let mut n = PospischilNeuron::new();
3223 n.step(f64::NAN);
3224 }
3225
3226 #[test]
3228 fn mainen_stable_without_input() {
3229 let mut n = MainenSejnowskiNeuron::new();
3232 for _ in 0..500 {
3233 n.step(0.0);
3234 }
3235 assert!(n.vs.is_finite());
3236 assert!(n.va.is_finite());
3237 }
3238 #[test]
3239 fn mainen_reset_clears_state() {
3240 let mut n = MainenSejnowskiNeuron::new();
3241 for _ in 0..100 {
3242 n.step(500.0);
3243 }
3244 n.reset();
3245 assert!((n.vs - (-65.0)).abs() < 1e-10);
3246 assert!((n.va - (-65.0)).abs() < 1e-10);
3247 }
3248 #[test]
3249 fn mainen_moderate_input_stable() {
3250 let mut n = MainenSejnowskiNeuron::new();
3252 for _ in 0..200 {
3253 n.step(500.0);
3254 }
3255 let _ = n.vs; }
3259 #[test]
3260 fn mainen_two_compartments_coupled() {
3261 let n = MainenSejnowskiNeuron::new();
3262 assert!(n.kappa > 0.0, "coupling should be positive");
3264 }
3265 #[test]
3266 fn mainen_weak_negative_no_crash() {
3267 let mut n = MainenSejnowskiNeuron::new();
3268 for _ in 0..200 {
3269 n.step(-10.0);
3270 }
3271 assert!(n.vs.is_finite());
3273 }
3274 #[test]
3275 fn mainen_nan_no_panic() {
3276 let mut n = MainenSejnowskiNeuron::new();
3277 n.step(f64::NAN);
3278 }
3279
3280 #[test]
3282 fn purkinje_silent_without_input() {
3283 let mut n = DeSchutterPurkinjeNeuron::new();
3284 let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3285 assert!(n.v.is_finite());
3287 }
3288 #[test]
3289 fn purkinje_reset_clears_state() {
3290 let mut n = DeSchutterPurkinjeNeuron::new();
3291 for _ in 0..100 {
3292 n.step(50.0);
3293 }
3294 n.reset();
3295 assert!((n.v - (-68.0)).abs() < 1e-10);
3296 assert!((n.ca - 0.0001).abs() < 1e-10);
3297 }
3298 #[test]
3299 fn purkinje_extreme_bounded() {
3300 let mut n = DeSchutterPurkinjeNeuron::new();
3301 for _ in 0..200 {
3302 n.step(1e4);
3303 }
3304 assert!(n.v.is_finite());
3305 }
3306 #[test]
3307 fn purkinje_ca_rises_with_spiking() {
3308 let mut n = DeSchutterPurkinjeNeuron::new();
3309 let ca_init = n.ca;
3310 for _ in 0..5000 {
3311 n.step(200.0);
3312 }
3313 assert!(n.ca > ca_init, "Ca²⁺ should rise during spiking");
3314 }
3315 #[test]
3316 fn purkinje_kca_activated_by_calcium() {
3317 let mut n = DeSchutterPurkinjeNeuron::new();
3318 for _ in 0..5000 {
3319 n.step(200.0);
3320 }
3321 assert!(
3322 n.q_kca > 0.0,
3323 "KCa should activate with Ca²⁺: q={}",
3324 n.q_kca
3325 );
3326 }
3327 #[test]
3328 fn purkinje_rk4_cross_backend_anchor() {
3329 let mut n = DeSchutterPurkinjeNeuron::new();
3330 let mut spikes = 0;
3331 for _ in 0..20_000 {
3332 spikes += n.step(500.0);
3333 }
3334 assert_eq!(spikes, 1);
3335 }
3336 #[test]
3337 fn purkinje_invalid_input_preserves_state() {
3338 let mut n = DeSchutterPurkinjeNeuron::new();
3339 for _ in 0..10 {
3340 let _ = n.step(200.0);
3341 }
3342 let old = (n.v, n.h_na, n.n_k, n.m_cap, n.h_cap, n.q_kca, n.ca);
3343 assert_eq!(n.step(f64::INFINITY), 0);
3344 assert_eq!((n.v, n.h_na, n.n_k, n.m_cap, n.h_cap, n.q_kca, n.ca), old);
3345 }
3346 #[test]
3347 fn purkinje_negative_no_crash() {
3348 let mut n = DeSchutterPurkinjeNeuron::new();
3349 for _ in 0..200 {
3350 n.step(-50.0);
3351 }
3352 assert!(n.v.is_finite());
3353 }
3354
3355 #[test]
3357 fn plant_r15_silent_without_input() {
3358 let mut n = PlantR15Neuron::new();
3359 for _ in 0..500 {
3361 n.step(0.0);
3362 }
3363 assert!(n.v.is_finite());
3364 }
3365 #[test]
3366 fn plant_r15_reset_clears_state() {
3367 let mut n = PlantR15Neuron::new();
3368 for _ in 0..100 {
3369 n.step(2.0);
3370 }
3371 n.reset();
3372 assert!((n.v - (-50.0)).abs() < 1e-10);
3373 assert!((n.ca - 0.1).abs() < 1e-10);
3374 }
3375 #[test]
3376 fn plant_r15_moderate_input_stable() {
3377 let mut n = PlantR15Neuron::new();
3379 for _ in 0..500 {
3380 n.step(2.0);
3381 }
3382 assert!(n.v.is_finite());
3383 }
3384 #[test]
3385 fn plant_r15_ca_dynamics() {
3386 let mut n = PlantR15Neuron::new();
3387 for _ in 0..500 {
3388 n.step(2.0);
3389 }
3390 assert!(n.ca >= 0.0, "Ca²⁺ must be non-negative");
3391 assert!(n.ca.is_finite());
3392 }
3393 #[test]
3394 fn plant_r15_weak_negative_no_crash() {
3395 let mut n = PlantR15Neuron::new();
3396 for _ in 0..200 {
3397 n.step(-1.0);
3398 }
3399 assert!(n.v.is_finite());
3400 }
3401 #[test]
3402 fn plant_r15_nan_no_panic() {
3403 let mut n = PlantR15Neuron::new();
3404 n.step(f64::NAN);
3405 }
3406
3407 #[test]
3409 fn prescott_zero_input_stable() {
3410 let mut n = PrescottNeuron::new();
3411 let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3412 assert!(n.v.is_finite());
3414 }
3415 #[test]
3416 fn prescott_reset_clears_state() {
3417 let mut n = PrescottNeuron::new();
3418 for _ in 0..100 {
3419 n.step(5.0);
3420 }
3421 n.reset();
3422 assert!((n.v - (-65.0)).abs() < 1e-10);
3423 assert!((n.w - 0.0).abs() < 1e-10);
3424 }
3425 #[test]
3426 fn prescott_rk4_reference_point() {
3427 let mut n = PrescottNeuron::new();
3428 assert_eq!(n.step(50.0), 0);
3429 assert!((n.v - (-44.498914201492525)).abs() < 1e-12);
3430 assert!((n.w - 1.4035864179018786e-05).abs() < 1e-17);
3431 }
3432 #[test]
3433 fn prescott_extreme_bounded() {
3434 let mut n = PrescottNeuron::new();
3435 for _ in 0..200 {
3436 n.step(1e4);
3437 }
3438 assert!(n.v.is_finite());
3439 }
3440 #[test]
3441 fn prescott_slow_var_adapts() {
3442 let mut n = PrescottNeuron::new();
3443 for _ in 0..500 {
3444 n.step(5.0);
3445 }
3446 assert!(n.w > 0.0, "slow variable w should activate during spiking");
3447 }
3448 #[test]
3449 fn prescott_negative_no_crash() {
3450 let mut n = PrescottNeuron::new();
3451 for _ in 0..200 {
3452 n.step(-10.0);
3453 }
3454 assert!(n.v.is_finite());
3455 }
3456 #[test]
3457 fn prescott_nan_no_panic() {
3458 let mut n = PrescottNeuron::new();
3459 n.step(f64::NAN);
3460 }
3461
3462 #[test]
3464 fn mn_silent_without_input() {
3465 let mut n = MihalasNieburNeuron::new();
3466 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
3467 assert_eq!(t, 0);
3468 }
3469 #[test]
3470 fn mn_reset_clears_state() {
3471 let mut n = MihalasNieburNeuron::new();
3472 for _ in 0..100 {
3473 n.step(5.0);
3474 }
3475 n.reset();
3476 assert!((n.v - n.v_rest).abs() < 1e-10);
3477 assert!((n.theta - n.theta_reset).abs() < 1e-10);
3478 }
3479 #[test]
3480 fn mn_rk4_reference_point() {
3481 let mut n = MihalasNieburNeuron::new();
3482 assert_eq!(n.step(0.5), 0);
3483 assert!((n.v - 0.04758125).abs() < 1e-12);
3484 assert!((n.theta - 1.0).abs() < 1e-15);
3485 assert!((n.i1 - 0.0).abs() < 1e-15);
3486 assert!((n.i2 - 0.0).abs() < 1e-15);
3487 }
3488 #[test]
3489 fn mn_spike_reset_uses_b() {
3490 let mut n = MihalasNieburNeuron::new();
3491 n.v = 0.99;
3492 n.b = 0.5;
3493 n.r1 = 1.25;
3494 n.r2 = -0.25;
3495 assert_eq!(n.step(2.0), 1);
3496 assert!((n.v - 0.5430570625).abs() < 1e-12);
3497 assert!((n.i1 - 1.25).abs() < 1e-15);
3498 assert!((n.i2 - (-0.25)).abs() < 1e-15);
3499 }
3500 #[test]
3501 fn mn_invalid_input_preserves_state() {
3502 let mut n = MihalasNieburNeuron::new();
3503 n.v = 0.2;
3504 n.i1 = 0.3;
3505 let before = (n.v, n.theta, n.i1, n.i2);
3506 assert_eq!(n.step(f64::NAN), 0);
3507 assert_eq!((n.v, n.theta, n.i1, n.i2), before);
3508 }
3509 #[test]
3510 fn mn_extreme_bounded() {
3511 let mut n = MihalasNieburNeuron::new();
3512 for _ in 0..200 {
3513 n.step(1e4);
3514 }
3515 assert!(n.v.is_finite());
3516 }
3517 #[test]
3518 fn mn_adaptive_threshold() {
3519 let mut n = MihalasNieburNeuron::new();
3520 n.a = 0.1;
3521 for _ in 0..100 {
3522 n.step(5.0);
3523 }
3524 assert!(n.theta.is_finite());
3526 }
3527 #[test]
3528 fn mn_negative_no_crash() {
3529 let mut n = MihalasNieburNeuron::new();
3530 for _ in 0..200 {
3531 n.step(-10.0);
3532 }
3533 assert!(n.v.is_finite());
3534 }
3535 #[test]
3536 fn mn_nan_no_panic() {
3537 let mut n = MihalasNieburNeuron::new();
3538 n.step(f64::NAN);
3539 }
3540
3541 #[test]
3543 fn glif_silent_without_input() {
3544 let mut n = GLIFNeuron::new();
3545 let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
3546 assert_eq!(t, 0);
3547 }
3548 #[test]
3549 fn glif_reset_clears_state() {
3550 let mut n = GLIFNeuron::new();
3551 for _ in 0..100 {
3552 n.step(30.0);
3553 }
3554 n.reset();
3555 assert!((n.v - n.v_rest).abs() < 1e-10);
3556 assert!((n.i_asc1).abs() < 1e-10);
3557 assert!((n.i_asc2).abs() < 1e-10);
3558 }
3559 #[test]
3560 fn glif_rk4_reference_point() {
3561 let mut n = GLIFNeuron::new();
3562 n.v = -68.0;
3563 n.theta = -45.0;
3564 n.i_asc1 = 0.4;
3565 n.i_asc2 = -0.2;
3566 assert_eq!(n.step(4.0), 0);
3567 assert!((n.v - (-67.7924658728125)).abs() < 1e-12);
3568 assert!((n.theta - (-45.049541282631253)).abs() < 1e-12);
3569 assert!((n.i_asc1 - 0.361935).abs() < 1e-12);
3570 assert!((n.i_asc2 - (-0.19900249583333334)).abs() < 1e-10);
3571 }
3572 #[test]
3573 fn glif_spike_reset_adds_candidate_threshold() {
3574 let mut n = GLIFNeuron::new();
3575 n.v = -51.0;
3576 n.theta = -50.5;
3577 n.delta_theta = 2.5;
3578 n.r_asc1 = 1.25;
3579 n.r_asc2 = -0.25;
3580 assert_eq!(n.step(40.0), 1);
3581 assert!((n.v - (-70.0)).abs() < 1e-12);
3582 assert!((n.theta - (-47.9930331381625)).abs() < 1e-12);
3583 assert!((n.i_asc1 - 1.25).abs() < 1e-12);
3584 assert!((n.i_asc2 - (-0.25)).abs() < 1e-12);
3585 }
3586 #[test]
3587 fn glif_invalid_input_preserves_state() {
3588 let mut n = GLIFNeuron::new();
3589 n.v = -68.0;
3590 n.i_asc1 = 0.4;
3591 let before = (n.v, n.theta, n.i_asc1, n.i_asc2);
3592 assert_eq!(n.step(f64::NAN), 0);
3593 assert_eq!((n.v, n.theta, n.i_asc1, n.i_asc2), before);
3594 }
3595 #[test]
3596 fn glif_extreme_bounded() {
3597 let mut n = GLIFNeuron::new();
3598 for _ in 0..200 {
3599 n.step(1e4);
3600 }
3601 assert!(n.v.is_finite());
3602 }
3603 #[test]
3604 fn glif_threshold_adapts_after_spike() {
3605 let mut n = GLIFNeuron::new();
3606 let theta_init = n.theta;
3607 for _ in 0..200 {
3608 n.step(30.0);
3609 }
3610 assert!(
3611 n.theta > theta_init,
3612 "theta should increase after spikes (delta_theta > 0)"
3613 );
3614 }
3615 #[test]
3616 fn glif_afterspike_currents() {
3617 let mut n = GLIFNeuron::new();
3618 for _ in 0..200 {
3619 n.step(30.0);
3620 }
3621 assert!(n.v.is_finite());
3623 }
3624 #[test]
3625 fn glif_negative_no_crash() {
3626 let mut n = GLIFNeuron::new();
3627 for _ in 0..200 {
3628 n.step(-30.0);
3629 }
3630 assert!(n.v.is_finite());
3631 }
3632
3633 #[test]
3635 fn gif_pop_exact_subthreshold_reference_point() {
3636 let mut n = GIFPopulationNeuron::new(7);
3637 n.v = -68.0;
3638 n.eta = 0.4;
3639 assert_eq!(n.step(4.0), 0);
3640 assert!((n.v - (-67.8370206677805)).abs() < 1e-12);
3641 assert!((n.eta - 0.398004991677073).abs() < 1e-15);
3642 }
3643 #[test]
3644 fn gif_pop_forced_spike_adds_decayed_adaptation() {
3645 let mut n = GIFPopulationNeuron::new(42);
3646 n.v = -51.0;
3647 n.eta = 0.3;
3648 n.theta = -90.0;
3649 n.lambda_0 = 1.0e9;
3650 assert_eq!(n.step(0.0), 1);
3651 assert!((n.v - n.v_reset).abs() < 1e-12);
3652 assert!((n.eta - 5.298503743757805).abs() < 1e-15);
3653 }
3654 #[test]
3655 fn gif_pop_invalid_input_preserves_state() {
3656 let mut n = GIFPopulationNeuron::new(42);
3657 n.v = -62.0;
3658 n.eta = 0.75;
3659 let before = (n.v, n.eta);
3660 assert_eq!(n.step(f64::NAN), 0);
3661 n.tau_m = 0.0;
3662 assert_eq!(n.step(1.0), 0);
3663 assert_eq!((n.v, n.eta), before);
3664 }
3665 #[test]
3666 fn gif_pop_seeded_reset_replays() {
3667 let mut n = GIFPopulationNeuron::new(123);
3668 n.theta = -90.0;
3669 n.lambda_0 = 1.0e9;
3670 let first: Vec<i32> = (0..3).map(|_| n.step(0.0)).collect();
3671 n.reset();
3672 let replay: Vec<i32> = (0..3).map(|_| n.step(0.0)).collect();
3673 assert_eq!(first, replay);
3674 assert!(n.eta > n.eta_increment);
3675 }
3676 #[test]
3677 fn gif_pop_negative_drive_remains_finite() {
3678 let mut n = GIFPopulationNeuron::new(42);
3679 for _ in 0..200 {
3680 n.step(-30.0);
3681 }
3682 assert!(n.v.is_finite());
3683 assert!(n.eta.is_finite());
3684 }
3685
3686 #[test]
3688 fn avron_rk4_reference_point() {
3689 let mut n = AvRonCardiacNeuron::new();
3690 n.v = -55.0;
3691 n.h = 0.55;
3692 n.n = 0.35;
3693 n.s = 0.45;
3694 assert_eq!(n.step(2.0), 0);
3695 assert!((n.v - (-50.0840498399381)).abs() < 1e-12);
3696 assert!((n.h - 0.5506609782132562).abs() < 1e-15);
3697 assert!((n.n - 0.34988677751350306).abs() < 1e-15);
3698 assert!((n.s - 0.4500091998827305).abs() < 1e-15);
3699 }
3700 #[test]
3701 fn avron_invalid_input_preserves_state() {
3702 let mut n = AvRonCardiacNeuron::new();
3703 n.v = -52.0;
3704 n.h = 0.4;
3705 n.n = 0.2;
3706 n.s = 0.8;
3707 let before = (n.v, n.h, n.n, n.s);
3708 assert_eq!(n.step(f64::NAN), 0);
3709 n.dt = 0.0;
3710 assert_eq!(n.step(1.0), 0);
3711 assert_eq!((n.v, n.h, n.n, n.s), before);
3712 }
3713 #[test]
3714 fn avron_rejects_corrupted_gate_state() {
3715 let mut n = AvRonCardiacNeuron::new();
3716 n.h = 1.2;
3717 let before = (n.v, n.h, n.n, n.s);
3718 assert_eq!(n.step(1.0), 0);
3719 assert_eq!((n.v, n.h, n.n, n.s), before);
3720 }
3721 #[test]
3722 fn avron_slow_current_remains_physical() {
3723 let mut n = AvRonCardiacNeuron::new();
3724 for _ in 0..1000 {
3725 n.step(5.0);
3726 assert!((0.0..=1.0).contains(&n.h));
3727 assert!((0.0..=1.0).contains(&n.n));
3728 assert!((0.0..=1.0).contains(&n.s));
3729 }
3730 assert!(n.v.is_finite());
3731 }
3732 #[test]
3733 fn avron_reset_clears_state() {
3734 let mut n = AvRonCardiacNeuron::new();
3735 n.v = -40.0;
3736 n.h = 0.1;
3737 n.n = 0.9;
3738 n.s = 0.2;
3739 n.reset();
3740 assert_eq!((n.v, n.h, n.n, n.s), (-60.0, 0.6, 0.3, 0.5));
3741 }
3742
3743 #[test]
3745 fn durstewitz_low_activity_zero_input() {
3746 let mut n = DurstewitzDopamineNeuron::new();
3747 let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3748 assert!(n.v.is_finite());
3750 }
3751 #[test]
3752 fn durstewitz_reset_clears_state() {
3753 let mut n = DurstewitzDopamineNeuron::new();
3754 for _ in 0..100 {
3755 n.step(3.0);
3756 }
3757 n.reset();
3758 assert!((n.v - (-65.0)).abs() < 1e-10);
3759 }
3760 #[test]
3761 fn durstewitz_extreme_bounded() {
3762 let mut n = DurstewitzDopamineNeuron::new();
3763 for _ in 0..200 {
3764 n.step(1e4);
3765 }
3766 assert!(n.v.is_finite());
3767 }
3768 #[test]
3769 fn durstewitz_d1_modulation() {
3770 let mut n_d1 = DurstewitzDopamineNeuron::new();
3772 n_d1.d1_level = 1.0;
3773 let mut n_no = DurstewitzDopamineNeuron::new();
3774 n_no.d1_level = 0.0;
3775 for _ in 0..1000 {
3776 n_d1.step(3.0);
3777 }
3778 for _ in 0..1000 {
3779 n_no.step(3.0);
3780 }
3781 assert!(n_d1.v.is_finite() && n_no.v.is_finite());
3783 }
3784 #[test]
3785 fn durstewitz_mg_block() {
3786 let n = DurstewitzDopamineNeuron::new();
3787 let block = 1.0 / (1.0 + n.mg * (-0.062 * n.v).exp() / 3.57);
3789 assert!(block < 0.1, "Mg²⁺ block at rest should be high: {}", block);
3790 }
3791 #[test]
3792 fn durstewitz_negative_no_crash() {
3793 let mut n = DurstewitzDopamineNeuron::new();
3794 for _ in 0..200 {
3795 n.step(-10.0);
3796 }
3797 assert!(n.v.is_finite());
3798 }
3799
3800 #[test]
3802 fn hill_tononi_silent_without_input() {
3803 let mut n = HillTononiNeuron::new();
3804 let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3805 assert!(n.v.is_finite());
3807 }
3808 #[test]
3809 fn hill_tononi_reset_clears_state() {
3810 let mut n = HillTononiNeuron::new();
3811 for _ in 0..100 {
3812 n.step(5.0);
3813 }
3814 n.reset();
3815 assert!((n.v - (-65.0)).abs() < 1e-10);
3816 assert!((n.na_i - 5.0).abs() < 1e-10);
3817 }
3818 #[test]
3819 fn hill_tononi_extreme_bounded() {
3820 let mut n = HillTononiNeuron::new();
3821 for _ in 0..200 {
3822 n.step(1e4);
3823 }
3824 assert!(n.v.is_finite());
3825 }
3826 #[test]
3827 fn hill_tononi_na_accumulation() {
3828 let mut n = HillTononiNeuron::new();
3829 for _ in 0..500 {
3830 n.step(5.0);
3831 }
3832 assert!(n.na_i.is_finite());
3834 assert!(n.na_i >= 0.0, "Na_i must be non-negative");
3835 }
3836 #[test]
3837 fn hill_tononi_negative_no_crash() {
3838 let mut n = HillTononiNeuron::new();
3839 for _ in 0..200 {
3840 n.step(-10.0);
3841 }
3842 assert!(n.v.is_finite());
3843 }
3844 #[test]
3845 fn hill_tononi_nan_no_panic() {
3846 let mut n = HillTononiNeuron::new();
3847 n.step(f64::NAN);
3848 }
3849
3850 #[test]
3852 fn bertram_silent_without_input() {
3853 let mut n = BertramPhantomBurster::new();
3854 for _ in 0..500 {
3855 n.step(0.0);
3856 }
3857 assert!(n.v.is_finite());
3858 }
3859 #[test]
3860 fn bertram_reset_clears_state() {
3861 let mut n = BertramPhantomBurster::new();
3862 for _ in 0..1000 {
3863 n.step(200.0);
3864 }
3865 n.reset();
3866 assert!((n.v - (-50.0)).abs() < 1e-10);
3867 assert!((n.s1 - 0.1).abs() < 1e-10);
3868 assert!((n.s2 - 0.1).abs() < 1e-10);
3869 }
3870 #[test]
3871 fn bertram_extreme_bounded() {
3872 let mut n = BertramPhantomBurster::new();
3873 for _ in 0..200 {
3874 n.step(1e4);
3875 }
3876 assert!(n.v.is_finite());
3877 }
3878 #[test]
3879 fn bertram_dual_slow_vars() {
3880 let mut n = BertramPhantomBurster::new();
3881 for _ in 0..10000 {
3882 n.step(200.0);
3883 }
3884 assert!(n.s1 >= 0.0 && n.s1 <= 1.0, "s1={}", n.s1);
3886 assert!(n.s2 >= 0.0 && n.s2 <= 1.0, "s2={}", n.s2);
3887 }
3888 #[test]
3889 fn bertram_negative_no_crash() {
3890 let mut n = BertramPhantomBurster::new();
3891 for _ in 0..200 {
3892 n.step(-100.0);
3893 }
3894 assert!(n.v.is_finite());
3895 }
3896 #[test]
3897 fn bertram_nan_no_panic() {
3898 let mut n = BertramPhantomBurster::new();
3899 n.step(f64::NAN);
3900 }
3901
3902 #[test]
3904 fn yamada_silent_without_input() {
3905 let mut n = YamadaNeuron::new();
3906 let t: i32 = (0..500).map(|_| n.step(0.0)).sum();
3907 assert_eq!(t, 0);
3908 }
3909 #[test]
3910 fn yamada_reset_clears_state() {
3911 let mut n = YamadaNeuron::new();
3912 for _ in 0..100 {
3913 n.step(5.0);
3914 }
3915 n.reset();
3916 assert!((n.v - (-60.0)).abs() < 1e-10);
3917 assert!((n.q - 0.0).abs() < 1e-10);
3918 }
3919 #[test]
3920 fn yamada_extreme_bounded() {
3921 let mut n = YamadaNeuron::new();
3922 for _ in 0..200 {
3923 n.step(1e4);
3924 }
3925 assert!(n.v.is_finite());
3926 }
3927 #[test]
3928 fn yamada_slow_q_adapts() {
3929 let mut n = YamadaNeuron::new();
3930 for _ in 0..2000 {
3931 n.step(5.0);
3932 }
3933 assert!(n.q > 0.0, "slow variable q should activate during spiking");
3934 }
3935 #[test]
3936 fn yamada_negative_no_crash() {
3937 let mut n = YamadaNeuron::new();
3938 for _ in 0..200 {
3939 n.step(-10.0);
3940 }
3941 assert!(n.v.is_finite());
3942 }
3943 #[test]
3944 fn yamada_nan_no_panic() {
3945 let mut n = YamadaNeuron::new();
3946 n.step(f64::NAN);
3947 }
3948}