sc_neurocore_engine 3.15.34

// SPDX-License-Identifier: AGPL-3.0-or-later
// Commercial license available
// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
// © Code 2020–2026 Miroslav Šotek. All rights reserved.
// ORCID: 0009-0009-3560-0851
// Contact: www.anulum.li | protoscience@anulum.li
// SC-NeuroCore — Hodgkin-Huxley-type conductance-based neuron models

//! Hodgkin-Huxley-type conductance-based neuron models.

use rand::{RngExt, SeedableRng};
use rand_xoshiro::Xoshiro256PlusPlus;

// ── Helper: safe alpha/beta kinetics avoiding division-by-zero ──

pub fn safe_rate(
    number_factor: f64,
    v_offset: f64,
    v: f64,
    denom_scale: f64,
    fallback: f64,
) -> f64 {
    let d = v + v_offset;
    if d.abs() < 1e-7 {
        fallback
    } else {
        number_factor * d / (1.0 - (-d / denom_scale).exp())
    }
}

/// Hodgkin-Huxley 1952 — 4-ODE ion channel model.
#[derive(Clone, Debug)]
pub struct HodgkinHuxleyNeuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub c_m: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl HodgkinHuxleyNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            m: 0.05,
            h: 0.6,
            n: 0.32,
            c_m: 1.0,
            g_na: 120.0,
            g_k: 36.0,
            g_l: 0.3,
            e_na: 50.0,
            e_k: -77.0,
            e_l: -54.4,
            dt: 0.01,
            v_threshold: 0.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let steps = (1.0 / self.dt) as usize;
        for _ in 0..steps {
            let am = safe_rate(0.1, 40.0, self.v, 10.0, 1.0);
            let bm = 4.0 * (-(self.v + 65.0) / 18.0).exp();
            let ah = 0.07 * (-(self.v + 65.0) / 20.0).exp();
            let bh = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
            let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
            let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -65.0;
        self.m = 0.05;
        self.h = 0.6;
        self.n = 0.32;
    }
}
impl Default for HodgkinHuxleyNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Traub-Miles — hippocampal CA3 pyramidal neuron with M-current.
///
/// Full Traub et al. 1991 model with Na, K_dr, leak, plus
/// M-current (Kv7/KCNQ) for spike frequency adaptation.
/// The M-current is the slow K⁺ conductance responsible for:
/// - Spike frequency adaptation (SFA)
/// - Theta-frequency resonance (~4-8 Hz)
/// - Subthreshold membrane potential oscillations
///
/// IM = g_M * w * (V - E_K)
/// dw/dt = (w_inf - w) / tau_w
/// w_inf = 1 / (1 + exp(-(V + 35) / 10))
/// tau_w = 100 / (3.3 * exp((V+35)/20) + exp(-(V+35)/20))
///
/// Traub et al., J Neurophysiol 66:635, 1991.
/// Yamada et al., Meth Neuronal Model (Koch & Segev), 1989 (M-current kinetics).
#[derive(Clone, Debug)]
pub struct TraubMilesNeuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl TraubMilesNeuron {
    pub fn new() -> Self {
        Self {
            v: -67.0,
            m: 0.05,
            h: 0.6,
            n: 0.3,
            g_na: 100.0,
            g_k: 80.0,
            g_l: 0.1,
            e_na: 50.0,
            e_k: -100.0,
            e_l: -67.0,
            dt: 0.01,
            v_threshold: -20.0,
        }
    }
    fn finite_gate(value: f64) -> bool {
        value.is_finite() && (0.0..=1.0).contains(&value)
    }
    fn valid_runtime(&self) -> bool {
        self.v.is_finite()
            && Self::finite_gate(self.m)
            && Self::finite_gate(self.h)
            && Self::finite_gate(self.n)
            && self.g_na.is_finite()
            && self.g_na >= 0.0
            && self.g_k.is_finite()
            && self.g_k >= 0.0
            && self.g_l.is_finite()
            && self.g_l >= 0.0
            && self.e_na.is_finite()
            && self.e_k.is_finite()
            && self.e_l.is_finite()
            && self.dt.is_finite()
            && self.dt > 0.0
            && self.v_threshold.is_finite()
    }
    fn rates(v: f64) -> Option<(f64, f64, f64, f64, f64, f64)> {
        let d = v + 54.0;
        let am = if d.abs() > 1e-6 {
            0.32 * d / (1.0 - (-d / 4.0).exp())
        } else {
            8.0
        };
        let d2 = v + 27.0;
        let bm = if d2.abs() > 1e-6 {
            0.28 * d2 / ((d2 / 5.0).exp() - 1.0)
        } else {
            5.6
        };
        let ah = 0.128 * (-(v + 50.0) / 18.0).exp();
        let bh = 4.0 / (1.0 + (-(v + 27.0) / 5.0).exp());
        let d3 = v + 52.0;
        let an = if d3.abs() > 1e-6 {
            0.032 * d3 / (1.0 - (-d3 / 5.0).exp())
        } else {
            0.32
        };
        let bn = 0.5 * (-(v + 57.0) / 40.0).exp();
        if [am, bm, ah, bh, an, bn]
            .iter()
            .all(|rate| rate.is_finite() && *rate >= 0.0)
        {
            Some((am, bm, ah, bh, an, bn))
        } else {
            None
        }
    }
    fn derivatives(
        &self,
        v: f64,
        m: f64,
        h: f64,
        n: f64,
        current: f64,
    ) -> Option<(f64, f64, f64, f64)> {
        if !v.is_finite() || !Self::finite_gate(m) || !Self::finite_gate(h) || !Self::finite_gate(n)
        {
            return None;
        }
        let (am, bm, ah, bh, an, bn) = Self::rates(v)?;
        let dm = am * (1.0 - m) - bm * m;
        let dh = ah * (1.0 - h) - bh * h;
        let dn = an * (1.0 - n) - bn * n;
        let i_na = self.g_na * m.powi(3) * h * (v - self.e_na);
        let i_k = self.g_k * n.powi(4) * (v - self.e_k);
        let i_l = self.g_l * (v - self.e_l);
        let dv = -i_na - i_k - i_l + current;
        if [dv, dm, dh, dn, i_na, i_k, i_l]
            .iter()
            .all(|value| value.is_finite())
        {
            Some((dv, dm, dh, dn))
        } else {
            None
        }
    }
    fn rk4_substep(
        &self,
        v: f64,
        m: f64,
        h: f64,
        n: f64,
        current: f64,
    ) -> Option<(f64, f64, f64, f64)> {
        let (k1_v, k1_m, k1_h, k1_n) = self.derivatives(v, m, h, n, current)?;
        let (k2_v, k2_m, k2_h, k2_n) = self.derivatives(
            v + 0.5 * self.dt * k1_v,
            m + 0.5 * self.dt * k1_m,
            h + 0.5 * self.dt * k1_h,
            n + 0.5 * self.dt * k1_n,
            current,
        )?;
        let (k3_v, k3_m, k3_h, k3_n) = self.derivatives(
            v + 0.5 * self.dt * k2_v,
            m + 0.5 * self.dt * k2_m,
            h + 0.5 * self.dt * k2_h,
            n + 0.5 * self.dt * k2_n,
            current,
        )?;
        let (k4_v, k4_m, k4_h, k4_n) = self.derivatives(
            v + self.dt * k3_v,
            m + self.dt * k3_m,
            h + self.dt * k3_h,
            n + self.dt * k3_n,
            current,
        )?;
        let next_v = v + self.dt * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v) / 6.0;
        let next_m = m + self.dt * (k1_m + 2.0 * k2_m + 2.0 * k3_m + k4_m) / 6.0;
        let next_h = h + self.dt * (k1_h + 2.0 * k2_h + 2.0 * k3_h + k4_h) / 6.0;
        let next_n = n + self.dt * (k1_n + 2.0 * k2_n + 2.0 * k3_n + k4_n) / 6.0;
        if next_v.is_finite()
            && Self::finite_gate(next_m)
            && Self::finite_gate(next_h)
            && Self::finite_gate(next_n)
        {
            Some((next_v, next_m, next_h, next_n))
        } else {
            None
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_runtime() {
            return 0;
        }
        let v_prev = self.v;
        let mut v = self.v;
        let mut m = self.m;
        let mut h = self.h;
        let mut n = self.n;
        for _ in 0..10 {
            let Some((next_v, next_m, next_h, next_n)) = self.rk4_substep(v, m, h, n, current)
            else {
                return 0;
            };
            v = next_v;
            m = next_m;
            h = next_h;
            n = next_n;
        }
        self.v = v;
        self.m = m;
        self.h = h;
        self.n = n;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        *self = Self::new();
    }
}
impl Default for TraubMilesNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Wang-Buzsaki — fast-spiking GABAergic interneuron. Wang & Buzsáki 1996.
#[derive(Clone, Debug)]
pub struct WangBuzsakiNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub phi: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl WangBuzsakiNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h: 0.8,
            n: 0.1,
            g_na: 35.0,
            g_k: 9.0,
            g_l: 0.1,
            e_na: 55.0,
            e_k: -90.0,
            e_l: -65.0,
            c_m: 1.0,
            phi: 5.0,
            dt: 0.01,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let n_sub = (0.5 / self.dt.max(0.001)) as usize;
        for _ in 0..n_sub {
            let am = safe_rate(0.1, 35.0, self.v, 10.0, 1.0);
            let bm = 4.0 * (-(self.v + 60.0) / 18.0).exp();
            let m_inf = am / (am + bm);
            let ah = 0.07 * (-(self.v + 58.0) / 20.0).exp();
            let bh = 1.0 / (1.0 + (-(self.v + 28.0) / 10.0).exp());
            let an = safe_rate(0.01, 34.0, self.v, 10.0, 0.1);
            let bn = 0.125 * (-(self.v + 44.0) / 80.0).exp();
            self.h += self.phi * (ah * (1.0 - self.h) - bh * self.h) * self.dt;
            self.n += self.phi * (an * (1.0 - self.n) - bn * self.n) * self.dt;
            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.v - self.e_na);
            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -65.0;
        self.h = 0.8;
        self.n = 0.1;
    }
}
impl Default for WangBuzsakiNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Connor-Stevens — A-type K current for delay tuning. Connor et al. 1977.
#[derive(Clone, Debug)]
pub struct ConnorStevensNeuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub a: f64,
    pub b: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_a: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_a: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl ConnorStevensNeuron {
    pub fn new() -> Self {
        Self {
            v: -68.0,
            m: 0.01,
            h: 0.99,
            n: 0.1,
            a: 0.5,
            b: 0.1,
            g_na: 120.0,
            g_k: 20.0,
            g_a: 47.7,
            g_l: 0.3,
            e_na: 55.0,
            e_k: -72.0,
            e_a: -75.0,
            e_l: -17.0,
            c_m: 1.0,
            dt: 0.01,
            v_threshold: 0.0,
        }
    }
    fn cs_safe_exp(x: f64) -> Option<f64> {
        if x.is_finite() && x <= 700.0 {
            Some(x.exp())
        } else {
            None
        }
    }

    fn cs_safe_rate(scale: f64, shift: f64, v: f64, denom: f64) -> Option<f64> {
        let delta = v + shift;
        let x = delta / denom;
        if x.abs() < 1e-9 {
            return Some(scale * denom);
        }
        let e = Self::cs_safe_exp(-x)?;
        let value = scale * delta / (1.0 - e);
        value.is_finite().then_some(value)
    }

    fn cs_valid_state(v: f64, m: f64, h: f64, n: f64, a: f64, b: f64) -> bool {
        [v, m, h, n, a, b].iter().all(|x| x.is_finite())
            && (-250.0..=250.0).contains(&v)
            && (-0.05..=1.05).contains(&m)
            && (-0.05..=1.05).contains(&h)
            && (-0.05..=1.05).contains(&n)
            && (-0.05..=1.5).contains(&a)
            && (-0.05..=1.05).contains(&b)
    }

    fn cs_valid_static(&self) -> bool {
        [
            self.g_na,
            self.g_k,
            self.g_a,
            self.g_l,
            self.e_na,
            self.e_k,
            self.e_a,
            self.e_l,
            self.c_m,
            self.dt,
            self.v_threshold,
        ]
        .iter()
        .all(|x| x.is_finite())
            && self.g_na >= 0.0
            && self.g_k >= 0.0
            && self.g_a >= 0.0
            && self.g_l >= 0.0
            && self.c_m > 0.0
            && self.dt > 0.0
    }

    fn cs_derivatives(
        &self,
        state: (f64, f64, f64, f64, f64, f64),
        current: f64,
    ) -> Option<(f64, f64, f64, f64, f64, f64)> {
        let (v, m, h, n, a, b) = state;
        let am = Self::cs_safe_rate(0.38, 29.7, v, 10.0)?;
        let bm = 15.2 * Self::cs_safe_exp(-(v + 54.7) / 18.0)?;
        let ah = 0.266 * Self::cs_safe_exp(-(v + 48.0) / 20.0)?;
        let bh = 3.8 / (1.0 + Self::cs_safe_exp(-(v + 18.0) / 10.0)?);
        let an = Self::cs_safe_rate(0.02, 45.7, v, 10.0)?;
        let bn = 0.25 * Self::cs_safe_exp(-(v + 55.7) / 80.0)?;
        let a_base = 0.0761 * Self::cs_safe_exp((v + 94.22) / 31.84)?
            / (1.0 + Self::cs_safe_exp((v + 1.17) / 28.93)?);
        if !a_base.is_finite() || a_base < 0.0 {
            return None;
        }
        let a_inf = a_base.powf(1.0 / 3.0);
        let tau_a = 0.3632 + 1.158 / (1.0 + Self::cs_safe_exp((v + 55.96) / 20.12)?);
        let b_base = 1.0 / (1.0 + Self::cs_safe_exp((v + 53.3) / 14.54)?);
        let b_inf = b_base.powf(4.0);
        let tau_b = 1.24 + 2.678 / (1.0 + Self::cs_safe_exp((v + 50.0) / 16.027)?);
        let i_na = self.g_na * m.powi(3) * h * (v - self.e_na);
        let i_k = self.g_k * n.powi(4) * (v - self.e_k);
        let i_a = self.g_a * a.powi(3) * b * (v - self.e_a);
        let i_l = self.g_l * (v - self.e_l);
        let deriv = (
            (-i_na - i_k - i_a - i_l + current) / self.c_m,
            am * (1.0 - m) - bm * m,
            ah * (1.0 - h) - bh * h,
            an * (1.0 - n) - bn * n,
            (a_inf - a) / tau_a,
            (b_inf - b) / tau_b,
        );
        [deriv.0, deriv.1, deriv.2, deriv.3, deriv.4, deriv.5]
            .iter()
            .all(|x| x.is_finite())
            .then_some(deriv)
    }

    fn cs_rk4_candidate(&self, current: f64) -> Option<(f64, f64, f64, f64, f64, f64)> {
        if !self.cs_valid_static()
            || !current.is_finite()
            || !Self::cs_valid_state(self.v, self.m, self.h, self.n, self.a, self.b)
        {
            return None;
        }
        let mut state = (self.v, self.m, self.h, self.n, self.a, self.b);
        let substeps = (1.0 / self.dt.max(0.001)) as usize;
        for _ in 0..substeps {
            let k1 = self.cs_derivatives(state, current)?;
            let k2_state = (
                state.0 + 0.5 * self.dt * k1.0,
                state.1 + 0.5 * self.dt * k1.1,
                state.2 + 0.5 * self.dt * k1.2,
                state.3 + 0.5 * self.dt * k1.3,
                state.4 + 0.5 * self.dt * k1.4,
                state.5 + 0.5 * self.dt * k1.5,
            );
            let k2 = self.cs_derivatives(k2_state, current)?;
            let k3_state = (
                state.0 + 0.5 * self.dt * k2.0,
                state.1 + 0.5 * self.dt * k2.1,
                state.2 + 0.5 * self.dt * k2.2,
                state.3 + 0.5 * self.dt * k2.3,
                state.4 + 0.5 * self.dt * k2.4,
                state.5 + 0.5 * self.dt * k2.5,
            );
            let k3 = self.cs_derivatives(k3_state, current)?;
            let k4_state = (
                state.0 + self.dt * k3.0,
                state.1 + self.dt * k3.1,
                state.2 + self.dt * k3.2,
                state.3 + self.dt * k3.3,
                state.4 + self.dt * k3.4,
                state.5 + self.dt * k3.5,
            );
            let k4 = self.cs_derivatives(k4_state, current)?;
            state = (
                state.0 + self.dt * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) / 6.0,
                state.1 + self.dt * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1) / 6.0,
                state.2 + self.dt * (k1.2 + 2.0 * k2.2 + 2.0 * k3.2 + k4.2) / 6.0,
                state.3 + self.dt * (k1.3 + 2.0 * k2.3 + 2.0 * k3.3 + k4.3) / 6.0,
                state.4 + self.dt * (k1.4 + 2.0 * k2.4 + 2.0 * k3.4 + k4.4) / 6.0,
                state.5 + self.dt * (k1.5 + 2.0 * k2.5 + 2.0 * k3.5 + k4.5) / 6.0,
            );
            if !Self::cs_valid_state(state.0, state.1, state.2, state.3, state.4, state.5) {
                return None;
            }
        }
        Some(state)
    }

    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let Some((v, m, h, n, a, b)) = self.cs_rk4_candidate(current) else {
            return 0;
        };
        self.v = v;
        self.m = m;
        self.h = h;
        self.n = n;
        self.a = a;
        self.b = b;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -68.0;
        self.m = 0.01;
        self.h = 0.99;
        self.n = 0.1;
        self.a = 0.5;
        self.b = 0.1;
    }
}
impl Default for ConnorStevensNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Destexhe thalamocortical relay neuron with T-current. Destexhe et al. 1993.
#[derive(Clone, Debug)]
pub struct DestexheThalamicNeuron {
    pub v: f64,
    pub h_na: f64,
    pub n_k: f64,
    pub m_t: f64,
    pub h_t: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_t: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_l: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl DestexheThalamicNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h_na: 0.6,
            n_k: 0.3,
            m_t: 0.0,
            h_t: 1.0,
            g_na: 100.0,
            g_k: 10.0,
            g_t: 2.0,
            g_l: 0.05,
            e_na: 50.0,
            e_k: -90.0,
            e_ca: 120.0,
            e_l: -70.0,
            dt: 0.02,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        for _ in 0..5 {
            let m_na = 1.0 / (1.0 + (-(self.v + 37.0) / 7.0).exp());
            let h_na_inf = 1.0 / (1.0 + ((self.v + 41.0) / 4.0).exp());
            let n_inf = 1.0 / (1.0 + (-(self.v + 25.0) / 12.0).exp());
            let m_t_inf = 1.0 / (1.0 + (-(self.v + 57.0) / 6.5).exp());
            let h_t_inf = 1.0 / (1.0 + ((self.v + 81.0) / 4.0).exp());
            // Voltage-dependent time constants (Destexhe 1993)
            let tau_h_na = (1.0
                / (0.128 * (-(self.v + 46.0) / 18.0).exp()
                    + 4.0 / (1.0 + (-(self.v + 23.0) / 5.0).exp())))
            .max(0.1);
            let tau_n_k = (1.0 / (0.032 * 5.0 + 0.5 * (-(self.v + 40.0) / 40.0).exp())).max(0.1);
            let tau_h_t = if self.v < -81.0 {
                (30.8
                    + 211.4 * ((self.v + 115.2) / 5.0).exp()
                        / (1.0 + ((self.v + 86.0) / 3.2).exp()))
                .max(0.1)
            } else {
                10.0
            };
            self.h_na += (h_na_inf - self.h_na) / tau_h_na * self.dt;
            self.n_k += (n_inf - self.n_k) / tau_n_k * self.dt;
            self.m_t = m_t_inf; // instantaneous (no ODE)
            self.h_t += (h_t_inf - self.h_t) / tau_h_t * self.dt;
            let i_na = self.g_na * m_na.powi(3) * self.h_na * (self.v - self.e_na);
            let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
            let i_t = self.g_t * self.m_t.powi(2) * self.h_t * (self.v - self.e_ca);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_t - i_l + current) * self.dt;
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -65.0;
        self.h_na = 0.6;
        self.n_k = 0.3;
        self.m_t = 0.0;
        self.h_t = 1.0;
    }
}
impl Default for DestexheThalamicNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Huber-Braun — temperature-sensitive cold receptor. Braun et al. 1998.
#[derive(Clone, Debug)]
pub struct HuberBraunNeuron {
    pub v: f64,
    pub a_sd: f64,
    pub a_sr: f64,
    pub g_sd: f64,
    pub g_sr: f64,
    pub g_l: f64,
    pub e_sd: f64,
    pub e_sr: f64,
    pub e_l: f64,
    pub tau_sd: f64,
    pub tau_sr: f64,
    pub eta: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl HuberBraunNeuron {
    pub fn new() -> Self {
        Self {
            v: -50.0,
            a_sd: 0.0,
            a_sr: 0.0,
            g_sd: 1.5,
            g_sr: 0.4,
            g_l: 0.1,
            e_sd: 50.0,
            e_sr: -90.0,
            e_l: -60.0,
            tau_sd: 10.0,
            tau_sr: 20.0,
            eta: 0.012,
            dt: 0.1,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        // Braun, Huber et al. 1998: SD V1/2 = -40 mV, slope = 6
        let sd_inf = 1.0 / (1.0 + (-(self.v + 40.0) / 6.0).exp());
        let sr_inf = 1.0 / (1.0 + ((self.v + 40.0) / 6.0).exp());
        self.a_sd += (sd_inf - self.a_sd) / self.tau_sd * self.dt;
        self.a_sr += (sr_inf - self.a_sr) / self.tau_sr * self.dt;
        let i_sd = self.g_sd * self.a_sd * (self.v - self.e_sd);
        let i_sr = self.g_sr * self.a_sr * (self.v - self.e_sr);
        let i_l = self.g_l * (self.v - self.e_l);
        self.v += (-i_sd - i_sr - i_l + current) * self.dt;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -50.0;
        self.a_sd = 0.0;
        self.a_sr = 0.0;
    }
}
impl Default for HuberBraunNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Golomb fast-spiking interneuron with Kv3. Golomb et al. 2007.
#[derive(Clone, Debug)]
pub struct GolombFSNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub p: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_kv3: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl GolombFSNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h: 0.9,
            n: 0.1,
            p: 0.0,
            g_na: 112.5,
            g_k: 225.0,
            g_kv3: 150.0,
            g_l: 0.25,
            e_na: 50.0,
            e_k: -90.0,
            e_l: -70.0,
            dt: 0.01,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        for _ in 0..10 {
            let m_inf = 1.0 / (1.0 + (-(self.v + 24.0) / 11.5).exp());
            let h_inf = 1.0 / (1.0 + ((self.v + 58.3) / 6.7).exp());
            let n_inf = 1.0 / (1.0 + (-(self.v + 12.4) / 6.8).exp());
            // Golomb et al. 2007: p_inf slope = 8.0
            let p_inf = 1.0 / (1.0 + (-(self.v + 3.0) / 8.0).exp());
            // Golomb et al. 2007: voltage-dependent time constants
            let tau_h = 0.5 + 14.0 / (1.0 + ((self.v + 60.0) / 12.0).exp());
            let tau_n = 0.087 + 11.4 / (1.0 + ((self.v + 14.6) / 8.6).exp());
            let tau_p = 0.1 + 4.0 / (1.0 + ((self.v + 25.0) / 10.0).exp());
            self.h += (h_inf - self.h) / tau_h * self.dt;
            self.n += (n_inf - self.n) / tau_n * self.dt;
            self.p += (p_inf - self.p) / tau_p * self.dt;
            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.v - self.e_na);
            // Golomb et al. 2007: I_Kd uses n^4, I_Kv3 uses p^2
            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
            let i_kv3 = self.g_kv3 * self.p.powi(2) * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_kv3 - i_l + current) * self.dt;
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -65.0;
        self.h = 0.9;
        self.n = 0.1;
        self.p = 0.0;
    }
}
impl Default for GolombFSNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Pospischil — minimal HH for 5 cortical cell types. Pospischil et al. 2008.
#[derive(Clone, Debug)]
pub struct PospischilNeuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub p: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_m: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub vt: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl PospischilNeuron {
    pub fn new() -> Self {
        Self {
            v: -70.0,
            m: 0.05,
            h: 0.6,
            n: 0.3,
            p: 0.0,
            g_na: 50.0,
            g_k: 5.0,
            g_m: 0.07,
            g_l: 0.1,
            e_na: 50.0,
            e_k: -90.0,
            e_l: -70.0,
            c_m: 1.0,
            vt: -56.2,
            dt: 0.025,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        for _ in 0..4 {
            let dv = self.v - self.vt;
            let x_m = dv - 13.0;
            let am = if x_m.abs() < 1e-6 {
                -0.32 * -4.0
            } else {
                -0.32 * x_m / ((-(x_m) / 4.0).exp() - 1.0)
            };
            let x_bm = dv - 40.0;
            let bm = if x_bm.abs() < 1e-6 {
                0.28 * 5.0
            } else {
                0.28 * x_bm / ((x_bm / 5.0).exp() - 1.0)
            };
            let ah = 0.128 * (-(dv - 17.0) / 18.0).exp();
            let bh = 4.0 / (1.0 + (-(dv - 40.0) / 5.0).exp());
            let x_n = dv - 15.0;
            let an = if x_n.abs() < 1e-6 {
                -0.032 * -5.0
            } else {
                -0.032 * x_n / ((-(x_n) / 5.0).exp() - 1.0)
            };
            let bn = 0.5 * (-(dv - 10.0) / 40.0).exp();
            let p_inf = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
            let tau_p =
                608.0 / (3.3 * ((self.v + 35.0) / 20.0).exp() + (-(self.v + 35.0) / 20.0).exp());
            self.p += (p_inf - self.p) / tau_p * self.dt;
            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
            let i_m = self.g_m * self.p * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_m - i_l + current) / self.c_m * self.dt;
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -70.0;
        self.m = 0.05;
        self.h = 0.6;
        self.n = 0.3;
        self.p = 0.0;
    }
}
impl Default for PospischilNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Mainen-Sejnowski — two-compartment (soma + axon). Mainen & Sejnowski 1996.
#[derive(Clone, Debug)]
pub struct MainenSejnowskiNeuron {
    pub vs: f64,
    pub va: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub kappa: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub c_s: f64,
    pub c_a: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl MainenSejnowskiNeuron {
    pub fn new() -> Self {
        Self {
            vs: -65.0,
            va: -65.0,
            m: 0.05,
            h: 0.6,
            n: 0.3,
            kappa: 10.0,
            g_na: 3000.0,
            g_k: 1500.0,
            g_l: 1.0,
            e_na: 50.0,
            e_k: -90.0,
            e_l: -70.0,
            c_s: 1.0,
            c_a: 0.1,
            dt: 0.005,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.vs;
        for _ in 0..20 {
            // Mainen & Sejnowski 1996 axonal rate functions
            let x_am = self.va + 25.0;
            let am = if x_am.abs() < 1e-6 {
                0.182 * 9.0
            } else {
                0.182 * x_am / (1.0 - (-(x_am) / 9.0).exp() + 1e-12)
            };
            let bm = if x_am.abs() < 1e-6 {
                0.124 * 9.0
            } else {
                -0.124 * x_am / (1.0 - ((x_am) / 9.0).exp() + 1e-12)
            };
            let x_ah = self.va + 40.0;
            let ah = if x_ah.abs() < 1e-6 {
                0.024 * 5.0
            } else {
                0.024 * x_ah / (1.0 - (-(x_ah) / 5.0).exp() + 1e-12)
            };
            let x_bh = self.va + 65.0;
            let bh = if x_bh.abs() < 1e-6 {
                0.0091 * 5.0
            } else {
                -0.0091 * x_bh / (1.0 - ((x_bh) / 5.0).exp() + 1e-12)
            };
            let x_an = self.va - 20.0;
            let an = if x_an.abs() < 1e-6 {
                0.02 * 9.0
            } else {
                0.02 * x_an / (1.0 - (-(x_an) / 9.0).exp() + 1e-12)
            };
            let bn = if x_an.abs() < 1e-6 {
                0.002 * 9.0
            } else {
                -0.002 * x_an / (1.0 - ((x_an) / 9.0).exp() + 1e-12)
            };
            self.m = (self.m + (am * (1.0 - self.m) - bm * self.m) * self.dt).clamp(0.0, 1.0);
            self.h = (self.h + (ah * (1.0 - self.h) - bh * self.h) * self.dt).clamp(0.0, 1.0);
            self.n = (self.n + (an * (1.0 - self.n) - bn * self.n) * self.dt).clamp(0.0, 1.0);
            let i_na = self.g_na * self.m.powi(3) * self.h * (self.va - self.e_na);
            let i_k = self.g_k * self.n * (self.va - self.e_k);
            let i_l_s = self.g_l * (self.vs - self.e_l);
            self.vs = (self.vs
                + (-i_l_s + self.kappa * (self.va - self.vs) + current) / self.c_s * self.dt)
                .clamp(-200.0, 200.0);
            self.va = (self.va
                + (-i_na - i_k + self.kappa * (self.vs - self.va)) / self.c_a * self.dt)
                .clamp(-200.0, 200.0);
        }
        if self.vs >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.vs = -65.0;
        self.va = -65.0;
        self.m = 0.05;
        self.h = 0.6;
        self.n = 0.3;
    }
}
impl Default for MainenSejnowskiNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// De Schutter-Bower Purkinje cell — Ca-dependent K. De Schutter & Bower 1994.
#[derive(Clone, Debug)]
pub struct DeSchutterPurkinjeNeuron {
    pub v: f64,
    pub h_na: f64,
    pub n_k: f64,
    pub m_cap: f64,
    pub h_cap: f64,
    pub q_kca: f64,
    pub ca: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_cap: f64,
    pub g_kca: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_l: f64,
    pub ca_decay: f64,
    pub f_ca: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl DeSchutterPurkinjeNeuron {
    pub fn new() -> Self {
        Self {
            v: -68.0,
            h_na: 0.8,
            n_k: 0.1,
            m_cap: 0.0,
            h_cap: 0.9,
            q_kca: 0.0,
            ca: 0.0001,
            g_na: 125.0,
            g_k: 10.0,
            g_cap: 45.0,
            g_kca: 35.0,
            g_l: 0.5,
            e_na: 45.0,
            e_k: -85.0,
            e_ca: 135.0,
            e_l: -68.0,
            ca_decay: 0.02,
            f_ca: 0.00024,
            dt: 0.01,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        for _ in 0..5 {
            let m_na = 1.0 / (1.0 + (-(self.v + 35.0) / 7.5).exp());
            let h_na_inf = 1.0 / (1.0 + ((self.v + 55.0) / 7.0).exp());
            let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 15.0).exp());
            let m_cap_inf = 1.0 / (1.0 + (-(self.v + 19.0) / 5.5).exp());
            let h_cap_inf = 1.0 / (1.0 + ((self.v + 48.0) / 7.0).exp());
            let q_inf = self.ca / (self.ca + 0.0002);
            let tau_h_na = 0.5 + 14.0 / (1.0 + ((self.v + 40.0) / 12.0).exp());
            let tau_n_k = 1.0 + 11.0 / (1.0 + ((self.v + 15.0) / 8.0).exp());
            self.h_na += (h_na_inf - self.h_na) / tau_h_na * self.dt;
            self.n_k += (n_inf - self.n_k) / tau_n_k * self.dt;
            self.m_cap += (m_cap_inf - self.m_cap) / 0.3 * self.dt;
            self.h_cap += (h_cap_inf - self.h_cap) / 45.0 * self.dt;
            self.q_kca += (q_inf - self.q_kca) / 1.0 * self.dt;
            let i_na = self.g_na * m_na.powi(3) * self.h_na * (self.v - self.e_na);
            let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
            let i_cap = self.g_cap * self.m_cap.powi(2) * self.h_cap * (self.v - self.e_ca);
            let i_kca = self.g_kca * self.q_kca * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_cap - i_kca - i_l + current) * self.dt;
            self.ca = (self.ca + (-self.f_ca * i_cap - self.ca_decay * self.ca) * self.dt).max(0.0);
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -68.0;
        self.h_na = 0.8;
        self.n_k = 0.1;
        self.m_cap = 0.0;
        self.h_cap = 0.9;
        self.q_kca = 0.0;
        self.ca = 0.0001;
    }
}
impl Default for DeSchutterPurkinjeNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Plant R15 — Aplysia parabolic burster. Plant & Kim 1976.
#[derive(Clone, Debug)]
pub struct PlantR15Neuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub ca: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_ca: f64,
    pub g_l: f64,
    pub g_kca: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub k_ca: f64,
    pub tau_ca: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl PlantR15Neuron {
    pub fn new() -> Self {
        Self {
            v: -50.0,
            m: 0.05,
            h: 0.6,
            n: 0.3,
            ca: 0.1,
            g_na: 4.0,
            g_k: 0.3,
            g_ca: 0.004,
            g_l: 0.003,
            g_kca: 0.03,
            e_na: 30.0,
            e_k: -75.0,
            e_ca: 140.0,
            e_l: -40.0,
            c_m: 1.0,
            k_ca: 0.0085,
            tau_ca: 500.0,
            dt: 0.05,
            v_threshold: -10.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        for _ in 0..5 {
            let am = safe_rate(0.1, 50.0, self.v, 10.0, 1.0);
            let bm = 4.0 * (-(self.v + 75.0) / 18.0).exp();
            let ah = 0.07 * (-(self.v + 50.0) / 20.0).exp();
            let bh = 1.0 / (1.0 + (-(self.v + 20.0) / 10.0).exp());
            let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
            let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
            let m_ca = 1.0 / (1.0 + (-(self.v + 25.0) / 5.0).exp());
            let kca_act = self.ca / (0.5 + self.ca);
            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
            let i_ca = self.g_ca * m_ca.powi(2) * (self.v - self.e_ca);
            let i_kca = self.g_kca * kca_act * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);
            self.v += (-i_na - i_k - i_ca - i_kca - i_l + current) / self.c_m * self.dt;
            self.ca = (self.ca + (-self.k_ca * i_ca - self.ca / self.tau_ca) * self.dt).max(0.0);
        }
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -50.0;
        self.m = 0.05;
        self.h = 0.6;
        self.n = 0.3;
        self.ca = 0.1;
    }
}
impl Default for PlantR15Neuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Prescott 2008 — Type I/II/III excitability tuning via M-current.
#[derive(Clone, Debug)]
pub struct PrescottNeuron {
    pub v: f64,
    pub w: f64,
    pub g_fast: f64,
    pub g_slow: f64,
    pub g_l: f64,
    pub e_fast: f64,
    pub e_slow: f64,
    pub e_l: f64,
    pub beta_w: f64,
    pub gamma_w: f64,
    pub tau_w: f64,
    pub phi: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl PrescottNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            w: 0.0,
            g_fast: 20.0,
            g_slow: 20.0,
            g_l: 2.0,
            e_fast: 50.0,
            e_slow: -100.0,
            e_l: -70.0,
            beta_w: -21.0,
            gamma_w: 15.0,
            tau_w: 100.0,
            phi: 0.15,
            dt: 0.1,
            v_threshold: -20.0,
        }
    }

    fn sigmoid(x: f64) -> f64 {
        if x >= 0.0 {
            let z = (-x).exp();
            1.0 / (1.0 + z)
        } else {
            let z = x.exp();
            z / (1.0 + z)
        }
    }

    fn valid_state(v: f64, w: f64) -> bool {
        v.is_finite() && w.is_finite() && (0.0..=1.0).contains(&w)
    }

    fn valid_runtime(&self) -> bool {
        Self::valid_state(self.v, self.w)
            && self.g_fast.is_finite()
            && self.g_fast >= 0.0
            && self.g_slow.is_finite()
            && self.g_slow >= 0.0
            && self.g_l.is_finite()
            && self.g_l >= 0.0
            && self.e_fast.is_finite()
            && self.e_slow.is_finite()
            && self.e_l.is_finite()
            && self.beta_w.is_finite()
            && self.gamma_w.is_finite()
            && self.gamma_w > 0.0
            && self.tau_w.is_finite()
            && self.tau_w > 0.0
            && self.phi.is_finite()
            && self.phi >= 0.0
            && self.dt.is_finite()
            && self.dt > 0.0
            && self.v_threshold.is_finite()
    }

    fn derivatives(&self, v: f64, w: f64, current: f64) -> Option<(f64, f64)> {
        if !Self::valid_state(v, w) {
            return None;
        }
        let m_inf = Self::sigmoid((v + 20.0) / 15.0);
        let w_inf = Self::sigmoid((v - self.beta_w) / self.gamma_w);
        let i_fast = self.g_fast * m_inf * (v - self.e_fast);
        let i_slow = self.g_slow * w * (v - self.e_slow);
        let i_l = self.g_l * (v - self.e_l);
        let dv = -i_fast - i_slow - i_l + current;
        let dw = self.phi * (w_inf - w) / self.tau_w;
        if dv.is_finite() && dw.is_finite() {
            Some((dv, dw))
        } else {
            None
        }
    }

    fn rk4_step(&self, current: f64) -> Option<(f64, f64)> {
        let dt = self.dt;
        let (k1_v, k1_w) = self.derivatives(self.v, self.w, current)?;
        let (k2_v, k2_w) =
            self.derivatives(self.v + 0.5 * dt * k1_v, self.w + 0.5 * dt * k1_w, current)?;
        let (k3_v, k3_w) =
            self.derivatives(self.v + 0.5 * dt * k2_v, self.w + 0.5 * dt * k2_w, current)?;
        let (k4_v, k4_w) = self.derivatives(self.v + dt * k3_v, self.w + dt * k3_w, current)?;
        let next_v = self.v + dt * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v) / 6.0;
        let next_w = self.w + dt * (k1_w + 2.0 * k2_w + 2.0 * k3_w + k4_w) / 6.0;
        if Self::valid_state(next_v, next_w) {
            Some((next_v, next_w))
        } else {
            None
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_runtime() {
            return 0;
        }
        let v_prev = self.v;
        let Some((next_v, next_w)) = self.rk4_step(current) else {
            return 0;
        };
        self.v = next_v;
        self.w = next_w;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }

    pub fn reset(&mut self) {
        self.v = -65.0;
        self.w = 0.0;
    }
}
impl Default for PrescottNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Mihalas-Niebur 2009 — generalized IF capturing 20 spike patterns.
#[derive(Clone, Debug)]
pub struct MihalasNieburNeuron {
    pub v: f64,
    pub theta: f64,
    pub i1: f64,
    pub i2: f64,
    pub v_rest: f64,
    pub v_reset: f64,
    pub theta_reset: f64,
    pub theta_inf: f64,
    pub tau_v: f64,
    pub tau_theta: f64,
    pub tau_1: f64,
    pub tau_2: f64,
    pub a: f64,
    pub b: f64,
    pub r1: f64,
    pub r2: f64,
    pub dt: f64,
}

impl MihalasNieburNeuron {
    pub fn new() -> Self {
        Self {
            v: 0.0,
            theta: 1.0,
            i1: 0.0,
            i2: 0.0,
            v_rest: 0.0,
            v_reset: 0.0,
            theta_reset: 1.0,
            theta_inf: 1.0,
            tau_v: 10.0,
            tau_theta: 100.0,
            tau_1: 10.0,
            tau_2: 200.0,
            a: 0.0,
            b: 0.0,
            r1: 0.0,
            r2: 0.0,
            dt: 1.0,
        }
    }

    fn finite_values(values: &[f64]) -> bool {
        values.iter().all(|value| value.is_finite())
    }

    fn valid_runtime(&self) -> bool {
        Self::finite_values(&[
            self.v,
            self.theta,
            self.i1,
            self.i2,
            self.v_rest,
            self.v_reset,
            self.theta_reset,
            self.theta_inf,
            self.tau_v,
            self.tau_theta,
            self.tau_1,
            self.tau_2,
            self.a,
            self.b,
            self.r1,
            self.r2,
            self.dt,
        ]) && self.tau_v > 0.0
            && self.tau_theta > 0.0
            && self.tau_1 > 0.0
            && self.tau_2 > 0.0
            && self.dt > 0.0
    }

    fn derivatives(&self, v: f64, theta: f64, i1: f64, i2: f64, current: f64) -> [f64; 4] {
        [
            (-(v - self.v_rest) + i1 + i2 + current) / self.tau_v,
            (self.theta_inf - theta + self.a * (v - self.v_rest)) / self.tau_theta,
            -i1 / self.tau_1,
            -i2 / self.tau_2,
        ]
    }

    fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
        [
            state[0] + scale * slope[0],
            state[1] + scale * slope[1],
            state[2] + scale * slope[2],
            state[3] + scale * slope[3],
        ]
    }

    fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
        let state = [self.v, self.theta, self.i1, self.i2];
        let half_dt = 0.5 * self.dt;
        let k1 = self.derivatives(state[0], state[1], state[2], state[3], current);
        let s2 = Self::add_scaled(state, k1, half_dt);
        let k2 = self.derivatives(s2[0], s2[1], s2[2], s2[3], current);
        let s3 = Self::add_scaled(state, k2, half_dt);
        let k3 = self.derivatives(s3[0], s3[1], s3[2], s3[3], current);
        let s4 = Self::add_scaled(state, k3, self.dt);
        let k4 = self.derivatives(s4[0], s4[1], s4[2], s4[3], current);
        let candidate = [
            state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
            state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
            state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
        ];
        if Self::finite_values(&candidate) {
            Some(candidate)
        } else {
            None
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_runtime() {
            return 0;
        }
        let Some(candidate) = self.rk4_candidate(current) else {
            return 0;
        };
        self.v = candidate[0];
        self.theta = candidate[1];
        self.i1 = candidate[2];
        self.i2 = candidate[3];
        if self.v >= self.theta {
            self.v = self.v_reset + self.b * (self.v - self.v_rest);
            self.theta = self.theta_reset.max(self.theta);
            self.i1 += self.r1;
            self.i2 += self.r2;
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = self.v_rest;
        self.theta = self.theta_reset;
        self.i1 = 0.0;
        self.i2 = 0.0;
    }
}
impl Default for MihalasNieburNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Allen GLIF5 — LIF + threshold adaptation + after-spike currents.
#[derive(Clone, Debug)]
pub struct GLIFNeuron {
    pub v: f64,
    pub theta: f64,
    pub theta_inf: f64,
    pub i_asc1: f64,
    pub i_asc2: f64,
    pub v_rest: f64,
    pub v_reset: f64,
    pub tau_m: f64,
    pub tau_theta: f64,
    pub tau_asc1: f64,
    pub tau_asc2: f64,
    pub a_theta: f64,
    pub delta_theta: f64,
    pub r_asc1: f64,
    pub r_asc2: f64,
    pub resistance: f64,
    pub dt: f64,
}

impl GLIFNeuron {
    pub fn new() -> Self {
        Self {
            v: -70.0,
            theta: -50.0,
            theta_inf: -50.0,
            i_asc1: 0.0,
            i_asc2: 0.0,
            v_rest: -70.0,
            v_reset: -70.0,
            tau_m: 10.0,
            tau_theta: 100.0,
            tau_asc1: 10.0,
            tau_asc2: 200.0,
            a_theta: 0.01,
            delta_theta: 2.0,
            r_asc1: 1.0,
            r_asc2: 0.5,
            resistance: 1.0,
            dt: 1.0,
        }
    }
    fn finite_values(values: &[f64]) -> bool {
        values.iter().all(|value| value.is_finite())
    }

    fn valid_runtime(&self) -> bool {
        Self::finite_values(&[
            self.v,
            self.theta,
            self.theta_inf,
            self.i_asc1,
            self.i_asc2,
            self.v_rest,
            self.v_reset,
            self.tau_m,
            self.tau_theta,
            self.tau_asc1,
            self.tau_asc2,
            self.a_theta,
            self.delta_theta,
            self.r_asc1,
            self.r_asc2,
            self.resistance,
            self.dt,
        ]) && self.tau_m > 0.0
            && self.tau_theta > 0.0
            && self.tau_asc1 > 0.0
            && self.tau_asc2 > 0.0
            && self.dt > 0.0
            && self.delta_theta >= 0.0
            && self.resistance >= 0.0
    }

    fn derivatives(&self, v: f64, theta: f64, i_asc1: f64, i_asc2: f64, current: f64) -> [f64; 4] {
        [
            (-(v - self.v_rest) + self.resistance * current + i_asc1 + i_asc2) / self.tau_m,
            (self.theta_inf - theta + self.a_theta * (v - self.v_rest)) / self.tau_theta,
            -i_asc1 / self.tau_asc1,
            -i_asc2 / self.tau_asc2,
        ]
    }

    fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
        [
            state[0] + scale * slope[0],
            state[1] + scale * slope[1],
            state[2] + scale * slope[2],
            state[3] + scale * slope[3],
        ]
    }

    fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
        let state = [self.v, self.theta, self.i_asc1, self.i_asc2];
        let half_dt = 0.5 * self.dt;
        let k1 = self.derivatives(state[0], state[1], state[2], state[3], current);
        let s2 = Self::add_scaled(state, k1, half_dt);
        let k2 = self.derivatives(s2[0], s2[1], s2[2], s2[3], current);
        let s3 = Self::add_scaled(state, k2, half_dt);
        let k3 = self.derivatives(s3[0], s3[1], s3[2], s3[3], current);
        let s4 = Self::add_scaled(state, k3, self.dt);
        let k4 = self.derivatives(s4[0], s4[1], s4[2], s4[3], current);
        let candidate = [
            state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
            state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
            state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
        ];
        if Self::finite_values(&candidate) {
            Some(candidate)
        } else {
            None
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_runtime() {
            return 0;
        }
        let Some(candidate) = self.rk4_candidate(current) else {
            return 0;
        };
        self.v = candidate[0];
        self.theta = candidate[1];
        self.i_asc1 = candidate[2];
        self.i_asc2 = candidate[3];
        if self.v >= self.theta {
            self.v = self.v_reset;
            self.theta += self.delta_theta;
            self.i_asc1 += self.r_asc1;
            self.i_asc2 += self.r_asc2;
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = self.v_rest;
        self.theta = self.theta_inf;
        self.i_asc1 = 0.0;
        self.i_asc2 = 0.0;
    }
}
impl Default for GLIFNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// GIF population — escape-rate generalized IF. Mensi et al. 2012.
#[derive(Clone, Debug)]
pub struct GIFPopulationNeuron {
    pub v: f64,
    pub theta: f64,
    pub eta: f64,
    pub tau_m: f64,
    pub tau_eta: f64,
    pub delta_v: f64,
    pub lambda_0: f64,
    pub eta_increment: f64,
    pub v_rest: f64,
    pub v_reset: f64,
    pub dt: f64,
    pub seed: u64,
    rng: Xoshiro256PlusPlus,
}

impl GIFPopulationNeuron {
    pub fn new(seed: u64) -> Self {
        Self {
            v: -65.0,
            theta: -50.0,
            eta: 0.0,
            tau_m: 20.0,
            tau_eta: 100.0,
            delta_v: 2.0,
            lambda_0: 0.001,
            eta_increment: 5.0,
            v_rest: -65.0,
            v_reset: -65.0,
            dt: 0.5,
            seed,
            rng: Xoshiro256PlusPlus::seed_from_u64(seed),
        }
    }

    fn finite_values(values: &[f64]) -> bool {
        values.iter().all(|value| value.is_finite())
    }

    fn valid_runtime(&self) -> bool {
        Self::finite_values(&[
            self.v,
            self.theta,
            self.eta,
            self.tau_m,
            self.tau_eta,
            self.delta_v,
            self.lambda_0,
            self.eta_increment,
            self.v_rest,
            self.v_reset,
            self.dt,
        ]) && self.tau_m > 0.0
            && self.tau_eta > 0.0
            && self.delta_v > 0.0
            && self.lambda_0 >= 0.0
            && self.dt > 0.0
    }

    fn advance_subthreshold(&self, current: f64) -> Option<(f64, f64)> {
        let eta_decay = (-self.dt / self.tau_eta).exp();
        let membrane_decay = (-self.dt / self.tau_m).exp();
        let x0 = self.v - self.v_rest - current;
        let eta_new = self.eta * eta_decay;
        let x_new = if (self.tau_m - self.tau_eta).abs() <= 1e-12 {
            membrane_decay * (x0 - self.eta * self.dt / self.tau_m)
        } else {
            let coupling = self.tau_eta / (self.tau_eta - self.tau_m);
            x0 * membrane_decay - self.eta * coupling * (eta_decay - membrane_decay)
        };
        let v_new = self.v_rest + current + x_new;
        if Self::finite_values(&[v_new, eta_new]) {
            Some((v_new, eta_new))
        } else {
            None
        }
    }

    fn spike_probability(&self, voltage: f64) -> f64 {
        if self.lambda_0 == 0.0 {
            return 0.0;
        }
        let exponent = ((voltage - self.theta) / self.delta_v).clamp(-745.0, 20.0);
        let hazard = self.lambda_0 * exponent.exp();
        (1.0 - (-hazard * self.dt).exp()).clamp(0.0, 1.0)
    }

    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_runtime() {
            return 0;
        }
        let Some((v_candidate, eta_candidate)) = self.advance_subthreshold(current) else {
            return 0;
        };
        self.v = v_candidate;
        self.eta = eta_candidate;
        if self.rng.random::<f64>() < self.spike_probability(self.v) {
            self.v = self.v_reset;
            self.eta += self.eta_increment;
            1
        } else {
            0
        }
    }

    pub fn reset(&mut self) {
        self.v = self.v_rest;
        self.eta = 0.0;
        self.rng = Xoshiro256PlusPlus::seed_from_u64(self.seed);
    }
}

/// Av-Ron cardiac ganglion — Type III burster. Av-Ron et al. 1991.
#[derive(Clone, Debug)]
pub struct AvRonCardiacNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub s: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_s: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_s: f64,
    pub e_l: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl AvRonCardiacNeuron {
    pub fn new() -> Self {
        Self {
            v: -60.0,
            h: 0.6,
            n: 0.3,
            s: 0.5,
            g_na: 80.0,
            g_k: 40.0,
            g_s: 20.0,
            g_l: 0.1,
            e_na: 40.0,
            e_k: -80.0,
            e_s: -25.0,
            e_l: -60.0,
            dt: 0.02,
            v_threshold: -20.0,
        }
    }

    fn finite_values(values: &[f64]) -> bool {
        values.iter().all(|value| value.is_finite())
    }

    fn gate_in_range(value: f64) -> bool {
        (0.0..=1.0).contains(&value)
    }

    fn bounded_exp(value: f64) -> f64 {
        value.clamp(-745.0, 709.0).exp()
    }

    fn sigmoid_pos(value: f64) -> f64 {
        1.0 / (1.0 + Self::bounded_exp(-value))
    }

    fn sigmoid_neg(value: f64) -> f64 {
        1.0 / (1.0 + Self::bounded_exp(value))
    }

    fn valid_runtime(&self) -> bool {
        Self::finite_values(&[
            self.v,
            self.h,
            self.n,
            self.s,
            self.g_na,
            self.g_k,
            self.g_s,
            self.g_l,
            self.e_na,
            self.e_k,
            self.e_s,
            self.e_l,
            self.dt,
            self.v_threshold,
        ]) && self.dt > 0.0
            && self.g_na >= 0.0
            && self.g_k >= 0.0
            && self.g_s >= 0.0
            && self.g_l >= 0.0
            && Self::gate_in_range(self.h)
            && Self::gate_in_range(self.n)
            && Self::gate_in_range(self.s)
    }

    fn rates(&self, voltage: f64) -> [f64; 7] {
        [
            Self::sigmoid_pos((voltage + 40.0) / 7.0),
            Self::sigmoid_neg((voltage + 45.0) / 5.0),
            Self::sigmoid_pos((voltage + 40.0) / 15.0),
            Self::sigmoid_neg((voltage + 35.0) / 3.0),
            1.0 + 12.0 * Self::sigmoid_neg((voltage + 50.0) / 8.0),
            1.0 + 8.0 * Self::sigmoid_neg((voltage + 35.0) / 8.0),
            200.0 + 1000.0 * Self::sigmoid_neg((voltage + 30.0) / 5.0),
        ]
    }

    fn derivatives(&self, state: [f64; 4], current: f64) -> [f64; 4] {
        let [voltage, h_gate, n_gate, s_gate] = state;
        if !Self::finite_values(&state)
            || !Self::gate_in_range(h_gate)
            || !Self::gate_in_range(n_gate)
            || !Self::gate_in_range(s_gate)
        {
            return [f64::NAN; 4];
        }
        let rates = self.rates(voltage);
        let i_na = self.g_na * rates[0].powi(3) * h_gate * (voltage - self.e_na);
        let i_k = self.g_k * n_gate.powi(4) * (voltage - self.e_k);
        let i_s = self.g_s * s_gate * (voltage - self.e_s);
        let i_l = self.g_l * (voltage - self.e_l);
        [
            -i_na - i_k - i_s - i_l + current,
            (rates[1] - h_gate) / rates[4],
            (rates[2] - n_gate) / rates[5],
            (rates[3] - s_gate) / rates[6],
        ]
    }

    fn add_scaled(state: [f64; 4], slope: [f64; 4], scale: f64) -> [f64; 4] {
        [
            state[0] + scale * slope[0],
            state[1] + scale * slope[1],
            state[2] + scale * slope[2],
            state[3] + scale * slope[3],
        ]
    }

    fn rk4_candidate(&self, current: f64) -> Option<[f64; 4]> {
        let state = [self.v, self.h, self.n, self.s];
        let half_dt = 0.5 * self.dt;
        let k1 = self.derivatives(state, current);
        let k2 = self.derivatives(Self::add_scaled(state, k1, half_dt), current);
        let k3 = self.derivatives(Self::add_scaled(state, k2, half_dt), current);
        let k4 = self.derivatives(Self::add_scaled(state, k3, self.dt), current);
        let candidate = [
            state[0] + self.dt * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            state[1] + self.dt * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0,
            state[2] + self.dt * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0,
            state[3] + self.dt * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]) / 6.0,
        ];
        if Self::finite_values(&candidate)
            && Self::gate_in_range(candidate[1])
            && Self::gate_in_range(candidate[2])
            && Self::gate_in_range(candidate[3])
        {
            Some(candidate)
        } else {
            None
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_runtime() {
            return 0;
        }
        let v_prev = self.v;
        let Some(candidate) = self.rk4_candidate(current) else {
            return 0;
        };
        self.v = candidate[0];
        self.h = candidate[1];
        self.n = candidate[2];
        self.s = candidate[3];
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }

    pub fn reset(&mut self) {
        self.v = -60.0;
        self.h = 0.6;
        self.n = 0.3;
        self.s = 0.5;
    }
}
impl Default for AvRonCardiacNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Durstewitz PFC neuron with D1 dopamine modulation. Durstewitz et al. 2000.
#[derive(Clone, Debug)]
pub struct DurstewitzDopamineNeuron {
    pub v: f64,
    pub h_na: f64,
    pub n_k: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_nmda: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_nmda: f64,
    pub e_l: f64,
    pub mg: f64,
    pub d1_level: f64,
    pub g_nmda_scale: f64,
    pub g_k_scale: f64,
    pub v_shift_na: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl DurstewitzDopamineNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h_na: 0.7,
            n_k: 0.2,
            g_na: 45.0,
            g_k: 18.0,
            g_nmda: 0.5,
            g_l: 0.02,
            e_na: 55.0,
            e_k: -80.0,
            e_nmda: 0.0,
            e_l: -65.0,
            mg: 1.0,
            d1_level: 0.0,
            g_nmda_scale: 2.5,
            g_k_scale: 1.5,
            v_shift_na: -5.0,
            dt: 0.05,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let shift = self.d1_level * self.v_shift_na;
        let m_inf = 1.0 / (1.0 + (-(self.v + 30.0 + shift) / 9.5).exp());
        let h_inf = 1.0 / (1.0 + ((self.v + 53.0) / 7.0).exp());
        let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 10.0).exp());
        let tau_h = 0.5 + 14.0 / (1.0 + ((self.v + 50.0) / 12.0).exp());
        let tau_n = 1.0 + 11.0 / (1.0 + ((self.v + 40.0) / 10.0).exp());
        self.h_na += (h_inf - self.h_na) / tau_h * self.dt;
        self.n_k += (n_inf - self.n_k) / tau_n * self.dt;
        let g_eff_nmda = self.g_nmda * (1.0 + self.d1_level * (self.g_nmda_scale - 1.0));
        let g_eff_k = self.g_k * (1.0 + self.d1_level * (self.g_k_scale - 1.0));
        let mg_block = 1.0 / (1.0 + self.mg * (-0.062 * self.v).exp() / 3.57);
        let i_na = self.g_na * m_inf.powi(3) * self.h_na * (self.v - self.e_na);
        let i_k = g_eff_k * self.n_k.powi(4) * (self.v - self.e_k);
        let i_nmda = g_eff_nmda * mg_block * (self.v - self.e_nmda);
        let i_l = self.g_l * (self.v - self.e_l);
        self.v += (-i_na - i_k - i_nmda - i_l + current) * self.dt;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -65.0;
        self.h_na = 0.7;
        self.n_k = 0.2;
    }
}
impl Default for DurstewitzDopamineNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Hill-Tononi 2005 — thalamocortical sleep/wake with Na-dependent K.
#[derive(Clone, Debug)]
pub struct HillTononiNeuron {
    pub v: f64,
    pub h_na: f64,
    pub n_k: f64,
    pub m_h: f64,
    pub h_t: f64,
    pub na_i: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_h: f64,
    pub g_t: f64,
    pub g_kna: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_h: f64,
    pub e_ca: f64,
    pub e_l: f64,
    pub na_pump_max: f64,
    pub na_eq: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl HillTononiNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h_na: 0.6,
            n_k: 0.3,
            m_h: 0.0,
            h_t: 0.9,
            na_i: 5.0,
            g_na: 50.0,
            g_k: 5.0,
            g_h: 1.0,
            g_t: 3.0,
            g_kna: 1.33,
            g_l: 0.02,
            e_na: 50.0,
            e_k: -90.0,
            e_h: -43.0,
            e_ca: 120.0,
            e_l: -70.0,
            na_pump_max: 20.0,
            na_eq: 9.5,
            dt: 0.05,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let m_na_inf = 1.0 / (1.0 + (-(self.v + 38.0) / 10.0).exp());
        let h_na_inf = 1.0 / (1.0 + ((self.v + 43.0) / 6.0).exp());
        let n_k_inf = 1.0 / (1.0 + (-(self.v + 27.0) / 11.5).exp());
        let m_h_inf = 1.0 / (1.0 + ((self.v + 75.0) / 5.5).exp());
        let m_t_inf = 1.0 / (1.0 + (-(self.v + 59.0) / 6.2).exp());
        let h_t_inf = 1.0 / (1.0 + ((self.v + 83.0) / 4.0).exp());
        let w_kna = 0.37 / (1.0 + (38.7 / self.na_i.max(0.01)).powf(3.5));
        let tau_h_na = (1.0 + 10.0 / (1.0 + ((self.v + 40.0) / 10.0).exp())).max(0.1);
        // Hill & Tononi 2005: tau_n_k = 5 + 47 * exp(-((V+50)/25)^2)
        let tau_n_k = (5.0 + 47.0 * (-(((self.v + 50.0) / 25.0).powi(2))).exp()).max(0.1);
        let tau_m_h = (20.0
            + 1000.0 / (((self.v + 71.5) / 14.2).exp() + (-(self.v + 89.0) / 11.6).exp()))
        .max(1.0);
        let tau_h_t = if self.v < -81.0 {
            (30.8 + 211.4 * ((self.v + 115.2) / 5.0).exp() / (1.0 + ((self.v + 86.0) / 3.2).exp()))
                .max(0.1)
        } else {
            10.0
        };
        self.h_na += (h_na_inf - self.h_na) / tau_h_na * self.dt;
        self.n_k += (n_k_inf - self.n_k) / tau_n_k * self.dt;
        self.m_h += (m_h_inf - self.m_h) / tau_m_h * self.dt;
        self.h_t += (h_t_inf - self.h_t) / tau_h_t * self.dt;
        let i_na = self.g_na * m_na_inf.powi(3) * self.h_na * (self.v - self.e_na);
        let i_k = self.g_k * self.n_k.powi(4) * (self.v - self.e_k);
        let i_h = self.g_h * self.m_h * (self.v - self.e_h);
        let i_t = self.g_t * m_t_inf.powi(2) * self.h_t * (self.v - self.e_ca);
        let i_kna = self.g_kna * w_kna * (self.v - self.e_k);
        let i_l = self.g_l * (self.v - self.e_l);
        self.v += (-i_na - i_k - i_h - i_t - i_kna - i_l + current) * self.dt;
        self.na_i = (self.na_i
            + (-0.001 * i_na - self.na_pump_max * (self.na_i / (self.na_i + self.na_eq)))
                * self.dt)
            .max(0.0);
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -65.0;
        self.h_na = 0.6;
        self.n_k = 0.3;
        self.m_h = 0.0;
        self.h_t = 0.9;
        self.na_i = 5.0;
    }
}
impl Default for HillTononiNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Bertram phantom burster — dual slow K for phantom bursting. Bertram et al. 2000.
#[derive(Clone, Debug)]
pub struct BertramPhantomBurster {
    pub v: f64,
    pub s1: f64,
    pub s2: f64,
    pub g_ca: f64,
    pub g_k: f64,
    pub g_s1: f64,
    pub g_s2: f64,
    pub g_l: f64,
    pub e_ca: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub v_m: f64,
    pub s_m: f64,
    pub v_n: f64,
    pub s_n: f64,
    pub v_s1: f64,
    pub s_s1: f64,
    pub v_s2: f64,
    pub s_s2: f64,
    pub tau_s1: f64,
    pub tau_s2: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl BertramPhantomBurster {
    pub fn new() -> Self {
        Self {
            v: -50.0,
            s1: 0.1,
            s2: 0.1,
            g_ca: 3.6,
            g_k: 10.0,
            g_s1: 4.0,
            g_s2: 4.0,
            g_l: 0.2,
            e_ca: 25.0,
            e_k: -75.0,
            e_l: -40.0,
            c_m: 5.3,
            v_m: -20.0,
            s_m: 12.0,
            v_n: -16.0,
            s_n: 5.6,
            v_s1: -40.0,
            s_s1: 10.0,
            v_s2: -42.0,
            s_s2: 0.4,
            tau_s1: 20000.0,
            tau_s2: 100000.0,
            dt: 0.5,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let m_inf = 1.0 / (1.0 + (-(self.v - self.v_m) / self.s_m).exp());
        let n_inf = 1.0 / (1.0 + (-(self.v - self.v_n) / self.s_n).exp());
        let s1_inf = 1.0 / (1.0 + (-(self.v - self.v_s1) / self.s_s1).exp());
        let s2_inf = 1.0 / (1.0 + (-(self.v - self.v_s2) / self.s_s2).exp());
        let i_ca = self.g_ca * m_inf * (self.v - self.e_ca);
        let i_k = self.g_k * n_inf * (self.v - self.e_k);
        let i_s1 = self.g_s1 * self.s1 * (self.v - self.e_k);
        let i_s2 = self.g_s2 * self.s2 * (self.v - self.e_k);
        let i_l = self.g_l * (self.v - self.e_l);
        self.v += (-i_ca - i_k - i_s1 - i_s2 - i_l + current) / self.c_m * self.dt;
        self.s1 += (s1_inf - self.s1) / self.tau_s1 * self.dt;
        self.s2 += (s2_inf - self.s2) / self.tau_s2 * self.dt;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -50.0;
        self.s1 = 0.1;
        self.s2 = 0.1;
    }
}
impl Default for BertramPhantomBurster {
    fn default() -> Self {
        Self::new()
    }
}

/// Yamada 1989 — subcritical Hopf burster.
#[derive(Clone, Debug)]
pub struct YamadaNeuron {
    pub v: f64,
    pub n: f64,
    pub q: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_q: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_q: f64,
    pub e_l: f64,
    pub tau_q: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl YamadaNeuron {
    pub fn new() -> Self {
        Self {
            v: -60.0,
            n: 0.1,
            q: 0.0,
            g_na: 20.0,
            g_k: 10.0,
            g_q: 5.0,
            g_l: 0.5,
            e_na: 60.0,
            e_k: -80.0,
            e_q: -80.0,
            e_l: -60.0,
            tau_q: 300.0,
            dt: 0.05,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let m_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 9.5).exp());
        let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 10.0).exp());
        let q_inf = 1.0 / (1.0 + (-(self.v + 50.0) / 10.0).exp());
        let tau_n = 1.0 + 7.5 / (1.0 + ((self.v + 40.0) / 12.0).exp());
        let i_na = self.g_na * m_inf.powi(3) * (1.0 - self.n) * (self.v - self.e_na);
        let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
        let i_q = self.g_q * self.q * (self.v - self.e_q);
        let i_l = self.g_l * (self.v - self.e_l);
        self.v += (-i_na - i_k - i_q - i_l + current) * self.dt;
        self.n += (n_inf - self.n) / tau_n * self.dt;
        self.q += (q_inf - self.q) / self.tau_q * self.dt;
        if self.v >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v = -60.0;
        self.n = 0.1;
        self.q = 0.0;
    }
}
impl Default for YamadaNeuron {
    fn default() -> Self {
        Self::new()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn hh_fires() {
        let mut n = HodgkinHuxleyNeuron::new();
        let t: i32 = (0..100).map(|_| n.step(10.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn traub_fires() {
        let mut n = TraubMilesNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn wb_fires() {
        let mut n = WangBuzsakiNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(2.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn cs_fires() {
        let mut n = ConnorStevensNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(10.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn destexhe_fires() {
        let mut n = DestexheThalamicNeuron::new();
        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn hb_fires() {
        let mut n = HuberBraunNeuron::new();
        let t: i32 = (0..5000).map(|_| n.step(10.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn golomb_fires() {
        let mut n = GolombFSNeuron::new();
        let t: i32 = (0..2000).map(|_| n.step(200.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn pospischil_fires() {
        let mut n = PospischilNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn mainen_fires() {
        let mut n = MainenSejnowskiNeuron::new();
        let t: i32 = (0..5000).map(|_| n.step(500.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn purkinje_fires() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        let t: i32 = (0..5000).map(|_| n.step(200.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn plant_r15_fires() {
        let mut n = PlantR15Neuron::new();
        let t: i32 = (0..500).map(|_| n.step(2.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn prescott_fires() {
        let mut n = PrescottNeuron::new();
        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn mn_fires() {
        let mut n = MihalasNieburNeuron::new();
        let t: i32 = (0..100).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn glif_fires() {
        let mut n = GLIFNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(30.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn gif_pop_fires() {
        let mut n = GIFPopulationNeuron::new(42);
        let t: i32 = (0..1000).map(|_| n.step(30.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn avron_fires() {
        let mut n = AvRonCardiacNeuron::new();
        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn durstewitz_fires() {
        let mut n = DurstewitzDopamineNeuron::new();
        let t: i32 = (0..1000).map(|_| n.step(3.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn hill_tononi_fires() {
        let mut n = HillTononiNeuron::new();
        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn bertram_fires() {
        let mut n = BertramPhantomBurster::new();
        let t: i32 = (0..10000).map(|_| n.step(200.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn yamada_fires() {
        let mut n = YamadaNeuron::new();
        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }

    // ── Multi-angle tests for all biophysical models ──

    // -- HodgkinHuxley --
    #[test]
    fn hh_silent_without_input() {
        let mut n = HodgkinHuxleyNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn hh_reset_clears_state() {
        let mut n = HodgkinHuxleyNeuron::new();
        for _ in 0..100 {
            n.step(10.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
        assert!((n.m - 0.05).abs() < 1e-10);
        assert!((n.h - 0.6).abs() < 1e-10);
        assert!((n.n - 0.32).abs() < 1e-10);
    }
    #[test]
    fn hh_extreme_input_bounded() {
        let mut n = HodgkinHuxleyNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn hh_gates_bounded() {
        let mut n = HodgkinHuxleyNeuron::new();
        for _ in 0..500 {
            n.step(10.0);
        }
        assert!(n.m >= 0.0 && n.m <= 1.0, "m={}", n.m);
        assert!(n.h >= 0.0 && n.h <= 1.0, "h={}", n.h);
        assert!(n.n >= 0.0 && n.n <= 1.0, "n={}", n.n);
    }
    #[test]
    fn hh_negative_input_no_crash() {
        let mut n = HodgkinHuxleyNeuron::new();
        for _ in 0..200 {
            n.step(-20.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn hh_nan_input_no_panic() {
        let mut n = HodgkinHuxleyNeuron::new();
        n.step(f64::NAN);
    }
    #[test]
    fn hh_sodium_potassium_opposition() {
        // Na activation drives depolarisation, K drives repolarisation
        let mut n = HodgkinHuxleyNeuron::new();
        for _ in 0..50 {
            n.step(10.0);
        }
        // After spiking, n (K activation) should have risen
        assert!(n.n > 0.32, "K activation n should increase during spiking");
    }

    // -- TraubMiles --
    #[test]
    fn traub_silent_without_input() {
        let mut n = TraubMilesNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn traub_reset_clears_state() {
        let mut n = TraubMilesNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-67.0)).abs() < 1e-10);
    }
    #[test]
    fn traub_extreme_bounded() {
        let mut n = TraubMilesNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn traub_gates_bounded() {
        let mut n = TraubMilesNeuron::new();
        for _ in 0..500 {
            n.step(5.0);
        }
        assert!(n.m >= 0.0 && n.m <= 1.01);
        assert!(n.h >= 0.0 && n.h <= 1.01);
        assert!(n.n >= 0.0 && n.n <= 1.01);
    }
    #[test]
    fn traub_weak_negative_no_crash() {
        let mut n = TraubMilesNeuron::new();
        for _ in 0..200 {
            n.step(-5.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn traub_nan_no_panic() {
        let mut n = TraubMilesNeuron::new();
        n.step(f64::NAN);
    }
    #[test]
    fn traub_rk4_reference_point() {
        let mut n = TraubMilesNeuron::new();
        n.v = -63.5;
        n.m = 0.08;
        n.h = 0.55;
        n.n = 0.32;
        let spike = n.step(4.0);
        assert_eq!(spike, 0);
        assert!((n.v - (-65.6638958700765)).abs() < 1e-13);
        assert!((n.m - 0.04237301812907925).abs() < 1e-15);
        assert!((n.h - 0.5626824931070477).abs() < 1e-15);
        assert!((n.n - 0.30356298261126924).abs() < 1e-15);
        assert!((n.v - (-65.66233161606698)).abs() > 1e-3);
    }

    // -- WangBuzsaki --
    #[test]
    fn wb_silent_without_input() {
        let mut n = WangBuzsakiNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn wb_reset_clears_state() {
        let mut n = WangBuzsakiNeuron::new();
        for _ in 0..100 {
            n.step(2.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn wb_extreme_bounded() {
        let mut n = WangBuzsakiNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn wb_fast_spiking_high_rate() {
        // WB model is fast-spiking — should achieve high rates
        let mut n = WangBuzsakiNeuron::new();
        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
        assert!(t > 10, "WB FS should produce many spikes, got {}", t);
    }
    #[test]
    fn wb_gates_bounded() {
        let mut n = WangBuzsakiNeuron::new();
        for _ in 0..500 {
            n.step(2.0);
        }
        assert!(n.h >= 0.0 && n.h <= 1.0);
        assert!(n.n >= 0.0 && n.n <= 1.0);
    }
    #[test]
    fn wb_negative_no_crash() {
        let mut n = WangBuzsakiNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn wb_nan_no_panic() {
        let mut n = WangBuzsakiNeuron::new();
        n.step(f64::NAN);
    }

    // -- ConnorStevens --
    #[test]
    fn cs_silent_without_input() {
        let mut n = ConnorStevensNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn cs_reset_clears_state() {
        let mut n = ConnorStevensNeuron::new();
        for _ in 0..100 {
            n.step(10.0);
        }
        n.reset();
        assert!((n.v - (-68.0)).abs() < 1e-10);
    }
    #[test]
    fn cs_extreme_bounded() {
        let mut n = ConnorStevensNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn cs_a_type_delays_spike() {
        // Connor-Stevens 1977: A-type K current causes onset delay
        // With 100 sub-steps/call, use more calls and verify A-type suppresses early firing
        let mut n_with_a = ConnorStevensNeuron::new();
        let mut n_no_a = ConnorStevensNeuron::new();
        n_no_a.g_a = 0.0;
        let spikes_with: i32 = (0..50).map(|_| n_with_a.step(8.0)).sum();
        let spikes_no: i32 = (0..50).map(|_| n_no_a.step(8.0)).sum();
        // Without A-current, neuron should fire more (no transient suppression)
        assert!(
            spikes_no >= spikes_with,
            "without A-type K, should fire more: no_a={} vs with_a={}",
            spikes_no,
            spikes_with
        );
    }
    #[test]
    fn cs_gates_bounded() {
        let mut n = ConnorStevensNeuron::new();
        for _ in 0..500 {
            n.step(10.0);
        }
        assert!(n.a >= 0.0 && n.a <= 1.5, "a={}", n.a); // a can slightly exceed 1 due to kinetics
        assert!(n.b >= 0.0 && n.b <= 1.0, "b={}", n.b);
    }
    #[test]
    fn cs_negative_no_crash() {
        let mut n = ConnorStevensNeuron::new();
        for _ in 0..200 {
            n.step(-20.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn cs_invalid_current_preserves_state() {
        let mut n = ConnorStevensNeuron::new();
        let before = (n.v, n.m, n.h, n.n, n.a, n.b);
        assert_eq!(n.step(f64::NAN), 0);
        assert_eq!((n.v, n.m, n.h, n.n, n.a, n.b), before);
    }
    #[test]
    fn cs_corrupt_runtime_state_preserves_state() {
        let mut n = ConnorStevensNeuron::new();
        n.b = f64::INFINITY;
        let before = (n.v, n.m, n.h, n.n, n.a, n.b);
        assert_eq!(n.step(6.0), 0);
        assert_eq!((n.v, n.m, n.h, n.n, n.a, n.b), before);
    }

    // -- Destexhe Thalamic --
    #[test]
    fn destexhe_no_crash_zero_input() {
        let mut n = DestexheThalamicNeuron::new();
        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        // Thalamic relays may have spontaneous activity via T-current
        assert!(n.v.is_finite());
    }
    #[test]
    fn destexhe_reset_clears_state() {
        let mut n = DestexheThalamicNeuron::new();
        for _ in 0..200 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn destexhe_extreme_bounded() {
        let mut n = DestexheThalamicNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn destexhe_t_current_rebound() {
        // T-type Ca²⁺ deinactivation: h_t_inf is high at hyperpolarised V
        // but voltage-dependent tau_h_t is very large (~37 s at V=-85),
        // so actual h_t recovery is slow — verify steady-state property instead.
        let v_hyp = -90.0_f64;
        let h_t_inf = 1.0 / (1.0 + ((v_hyp + 81.0) / 4.0).exp());
        assert!(h_t_inf > 0.9, "h_t_inf should be ~1 at V=-90: {}", h_t_inf);

        // Verify model stability through hyperpolarise-release cycle
        let mut n = DestexheThalamicNeuron::new();
        for _ in 0..500 {
            n.step(-5.0);
        }
        let v_before_release = n.v;
        for _ in 0..500 {
            n.step(0.0);
        }
        assert!(n.v.is_finite());
        assert!(
            n.v > v_before_release,
            "V should increase after release (rebound)"
        );
    }
    #[test]
    fn destexhe_negative_no_crash() {
        let mut n = DestexheThalamicNeuron::new();
        for _ in 0..200 {
            n.step(-20.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn destexhe_nan_no_panic() {
        let mut n = DestexheThalamicNeuron::new();
        n.step(f64::NAN);
    }

    // -- HuberBraun --
    #[test]
    fn hb_silent_without_input() {
        let mut n = HuberBraunNeuron::new();
        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        // HuberBraun may have spontaneous activity, so just check finite
        assert!(n.v.is_finite());
    }
    #[test]
    fn hb_reset_clears_state() {
        let mut n = HuberBraunNeuron::new();
        for _ in 0..200 {
            n.step(10.0);
        }
        n.reset();
        assert!((n.v - (-50.0)).abs() < 1e-10);
    }
    #[test]
    fn hb_extreme_bounded() {
        let mut n = HuberBraunNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn hb_negative_no_crash() {
        let mut n = HuberBraunNeuron::new();
        for _ in 0..500 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn hb_nan_no_panic() {
        let mut n = HuberBraunNeuron::new();
        n.step(f64::NAN);
    }

    // -- GolombFS --
    #[test]
    fn golomb_silent_without_input() {
        let mut n = GolombFSNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn golomb_reset_clears_state() {
        let mut n = GolombFSNeuron::new();
        for _ in 0..100 {
            n.step(200.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn golomb_extreme_bounded() {
        // Golomb et al. 2007: n^4 kinetics diverge at extreme I; test at high but realistic drive
        let mut n = GolombFSNeuron::new();
        for _ in 0..200 {
            n.step(200.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn golomb_kv3_enables_fast_spiking() {
        // Kv3 current enables high-frequency firing
        let mut n = GolombFSNeuron::new();
        let t: i32 = (0..5000).map(|_| n.step(300.0)).sum();
        assert!(t > 0, "Golomb FS should fire with high input, got {}", t);
    }
    #[test]
    fn golomb_negative_no_crash() {
        let mut n = GolombFSNeuron::new();
        for _ in 0..200 {
            n.step(-100.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn golomb_nan_no_panic() {
        let mut n = GolombFSNeuron::new();
        n.step(f64::NAN);
    }

    // -- Pospischil --
    #[test]
    fn pospischil_silent_without_input() {
        let mut n = PospischilNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn pospischil_reset_clears_state() {
        let mut n = PospischilNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-70.0)).abs() < 1e-10);
    }
    #[test]
    fn pospischil_moderate_input_stable() {
        let mut n = PospischilNeuron::new();
        for _ in 0..200 {
            n.step(10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn pospischil_m_current_present() {
        let mut n = PospischilNeuron::new();
        for _ in 0..200 {
            n.step(5.0);
        }
        assert!(n.p > 0.0, "M-current (p) should activate during spiking");
    }
    #[test]
    fn pospischil_negative_no_crash() {
        let mut n = PospischilNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn pospischil_nan_no_panic() {
        let mut n = PospischilNeuron::new();
        n.step(f64::NAN);
    }

    // -- MainenSejnowski --
    #[test]
    fn mainen_stable_without_input() {
        // Mainen 1996 model may produce transient spikes at I=0
        // (confirmed in Python reference). Verify stability only.
        let mut n = MainenSejnowskiNeuron::new();
        for _ in 0..500 {
            n.step(0.0);
        }
        assert!(n.vs.is_finite());
        assert!(n.va.is_finite());
    }
    #[test]
    fn mainen_reset_clears_state() {
        let mut n = MainenSejnowskiNeuron::new();
        for _ in 0..100 {
            n.step(500.0);
        }
        n.reset();
        assert!((n.vs - (-65.0)).abs() < 1e-10);
        assert!((n.va - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn mainen_moderate_input_stable() {
        // Two-compartment model with high conductances — moderate input
        let mut n = MainenSejnowskiNeuron::new();
        for _ in 0..200 {
            n.step(500.0);
        }
        // High-conductance 2-compartment may diverge at extremes;
        // test moderate stability
        let _ = n.vs; // no panic
    }
    #[test]
    fn mainen_two_compartments_coupled() {
        let n = MainenSejnowskiNeuron::new();
        // kappa > 0 means compartments are coupled
        assert!(n.kappa > 0.0, "coupling should be positive");
    }
    #[test]
    fn mainen_weak_negative_no_crash() {
        let mut n = MainenSejnowskiNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        // Weak negative is safer for 2-compartment
        assert!(n.vs.is_finite());
    }
    #[test]
    fn mainen_nan_no_panic() {
        let mut n = MainenSejnowskiNeuron::new();
        n.step(f64::NAN);
    }

    // -- DeSchutterPurkinje --
    #[test]
    fn purkinje_silent_without_input() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        // Purkinje cells may have some spontaneous activity
        assert!(n.v.is_finite());
    }
    #[test]
    fn purkinje_reset_clears_state() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        for _ in 0..100 {
            n.step(50.0);
        }
        n.reset();
        assert!((n.v - (-68.0)).abs() < 1e-10);
        assert!((n.ca - 0.0001).abs() < 1e-10);
    }
    #[test]
    fn purkinje_extreme_bounded() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn purkinje_ca_rises_with_spiking() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        let ca_init = n.ca;
        for _ in 0..5000 {
            n.step(200.0);
        }
        assert!(n.ca > ca_init, "Ca²⁺ should rise during spiking");
    }
    #[test]
    fn purkinje_kca_activated_by_calcium() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        for _ in 0..5000 {
            n.step(200.0);
        }
        assert!(
            n.q_kca > 0.0,
            "KCa should activate with Ca²⁺: q={}",
            n.q_kca
        );
    }
    #[test]
    fn purkinje_negative_no_crash() {
        let mut n = DeSchutterPurkinjeNeuron::new();
        for _ in 0..200 {
            n.step(-50.0);
        }
        assert!(n.v.is_finite());
    }

    // -- PlantR15 --
    #[test]
    fn plant_r15_silent_without_input() {
        let mut n = PlantR15Neuron::new();
        // R15 is a parabolic burster — may burst spontaneously
        for _ in 0..500 {
            n.step(0.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn plant_r15_reset_clears_state() {
        let mut n = PlantR15Neuron::new();
        for _ in 0..100 {
            n.step(2.0);
        }
        n.reset();
        assert!((n.v - (-50.0)).abs() < 1e-10);
        assert!((n.ca - 0.1).abs() < 1e-10);
    }
    #[test]
    fn plant_r15_moderate_input_stable() {
        // Plant R15 is a parabolic burster — moderate input stability
        let mut n = PlantR15Neuron::new();
        for _ in 0..500 {
            n.step(2.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn plant_r15_ca_dynamics() {
        let mut n = PlantR15Neuron::new();
        for _ in 0..500 {
            n.step(2.0);
        }
        assert!(n.ca >= 0.0, "Ca²⁺ must be non-negative");
        assert!(n.ca.is_finite());
    }
    #[test]
    fn plant_r15_weak_negative_no_crash() {
        let mut n = PlantR15Neuron::new();
        for _ in 0..200 {
            n.step(-1.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn plant_r15_nan_no_panic() {
        let mut n = PlantR15Neuron::new();
        n.step(f64::NAN);
    }

    // -- Prescott --
    #[test]
    fn prescott_zero_input_stable() {
        let mut n = PrescottNeuron::new();
        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        // Prescott has fast Na conductance — may produce spontaneous activity
        assert!(n.v.is_finite());
    }
    #[test]
    fn prescott_reset_clears_state() {
        let mut n = PrescottNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
        assert!((n.w - 0.0).abs() < 1e-10);
    }
    #[test]
    fn prescott_rk4_reference_point() {
        let mut n = PrescottNeuron::new();
        assert_eq!(n.step(50.0), 0);
        assert!((n.v - (-44.498914201492525)).abs() < 1e-12);
        assert!((n.w - 1.4035864179018786e-05).abs() < 1e-17);
    }
    #[test]
    fn prescott_extreme_bounded() {
        let mut n = PrescottNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn prescott_slow_var_adapts() {
        let mut n = PrescottNeuron::new();
        for _ in 0..500 {
            n.step(5.0);
        }
        assert!(n.w > 0.0, "slow variable w should activate during spiking");
    }
    #[test]
    fn prescott_negative_no_crash() {
        let mut n = PrescottNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn prescott_nan_no_panic() {
        let mut n = PrescottNeuron::new();
        n.step(f64::NAN);
    }

    // -- MihalasNiebur --
    #[test]
    fn mn_silent_without_input() {
        let mut n = MihalasNieburNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn mn_reset_clears_state() {
        let mut n = MihalasNieburNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - n.v_rest).abs() < 1e-10);
        assert!((n.theta - n.theta_reset).abs() < 1e-10);
    }
    #[test]
    fn mn_rk4_reference_point() {
        let mut n = MihalasNieburNeuron::new();
        assert_eq!(n.step(0.5), 0);
        assert!((n.v - 0.04758125).abs() < 1e-12);
        assert!((n.theta - 1.0).abs() < 1e-15);
        assert!((n.i1 - 0.0).abs() < 1e-15);
        assert!((n.i2 - 0.0).abs() < 1e-15);
    }
    #[test]
    fn mn_spike_reset_uses_b() {
        let mut n = MihalasNieburNeuron::new();
        n.v = 0.99;
        n.b = 0.5;
        n.r1 = 1.25;
        n.r2 = -0.25;
        assert_eq!(n.step(2.0), 1);
        assert!((n.v - 0.5430570625).abs() < 1e-12);
        assert!((n.i1 - 1.25).abs() < 1e-15);
        assert!((n.i2 - (-0.25)).abs() < 1e-15);
    }
    #[test]
    fn mn_invalid_input_preserves_state() {
        let mut n = MihalasNieburNeuron::new();
        n.v = 0.2;
        n.i1 = 0.3;
        let before = (n.v, n.theta, n.i1, n.i2);
        assert_eq!(n.step(f64::NAN), 0);
        assert_eq!((n.v, n.theta, n.i1, n.i2), before);
    }
    #[test]
    fn mn_extreme_bounded() {
        let mut n = MihalasNieburNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn mn_adaptive_threshold() {
        let mut n = MihalasNieburNeuron::new();
        n.a = 0.1;
        for _ in 0..100 {
            n.step(5.0);
        }
        // Threshold should have adapted
        assert!(n.theta.is_finite());
    }
    #[test]
    fn mn_negative_no_crash() {
        let mut n = MihalasNieburNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn mn_nan_no_panic() {
        let mut n = MihalasNieburNeuron::new();
        n.step(f64::NAN);
    }

    // -- GLIF --
    #[test]
    fn glif_silent_without_input() {
        let mut n = GLIFNeuron::new();
        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn glif_reset_clears_state() {
        let mut n = GLIFNeuron::new();
        for _ in 0..100 {
            n.step(30.0);
        }
        n.reset();
        assert!((n.v - n.v_rest).abs() < 1e-10);
        assert!((n.i_asc1).abs() < 1e-10);
        assert!((n.i_asc2).abs() < 1e-10);
    }
    #[test]
    fn glif_rk4_reference_point() {
        let mut n = GLIFNeuron::new();
        n.v = -68.0;
        n.theta = -45.0;
        n.i_asc1 = 0.4;
        n.i_asc2 = -0.2;
        assert_eq!(n.step(4.0), 0);
        assert!((n.v - (-67.7924658728125)).abs() < 1e-12);
        assert!((n.theta - (-45.049541282631253)).abs() < 1e-12);
        assert!((n.i_asc1 - 0.361935).abs() < 1e-12);
        assert!((n.i_asc2 - (-0.19900249583333334)).abs() < 1e-10);
    }
    #[test]
    fn glif_spike_reset_adds_candidate_threshold() {
        let mut n = GLIFNeuron::new();
        n.v = -51.0;
        n.theta = -50.5;
        n.delta_theta = 2.5;
        n.r_asc1 = 1.25;
        n.r_asc2 = -0.25;
        assert_eq!(n.step(40.0), 1);
        assert!((n.v - (-70.0)).abs() < 1e-12);
        assert!((n.theta - (-47.9930331381625)).abs() < 1e-12);
        assert!((n.i_asc1 - 1.25).abs() < 1e-12);
        assert!((n.i_asc2 - (-0.25)).abs() < 1e-12);
    }
    #[test]
    fn glif_invalid_input_preserves_state() {
        let mut n = GLIFNeuron::new();
        n.v = -68.0;
        n.i_asc1 = 0.4;
        let before = (n.v, n.theta, n.i_asc1, n.i_asc2);
        assert_eq!(n.step(f64::NAN), 0);
        assert_eq!((n.v, n.theta, n.i_asc1, n.i_asc2), before);
    }
    #[test]
    fn glif_extreme_bounded() {
        let mut n = GLIFNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn glif_threshold_adapts_after_spike() {
        let mut n = GLIFNeuron::new();
        let theta_init = n.theta;
        for _ in 0..200 {
            n.step(30.0);
        }
        assert!(
            n.theta > theta_init,
            "theta should increase after spikes (delta_theta > 0)"
        );
    }
    #[test]
    fn glif_afterspike_currents() {
        let mut n = GLIFNeuron::new();
        for _ in 0..200 {
            n.step(30.0);
        }
        // After spiking, ASC should have been triggered (then decayed)
        assert!(n.v.is_finite());
    }
    #[test]
    fn glif_negative_no_crash() {
        let mut n = GLIFNeuron::new();
        for _ in 0..200 {
            n.step(-30.0);
        }
        assert!(n.v.is_finite());
    }

    // -- GIFPopulation --
    #[test]
    fn gif_pop_exact_subthreshold_reference_point() {
        let mut n = GIFPopulationNeuron::new(7);
        n.v = -68.0;
        n.eta = 0.4;
        assert_eq!(n.step(4.0), 0);
        assert!((n.v - (-67.8370206677805)).abs() < 1e-12);
        assert!((n.eta - 0.398004991677073).abs() < 1e-15);
    }
    #[test]
    fn gif_pop_forced_spike_adds_decayed_adaptation() {
        let mut n = GIFPopulationNeuron::new(42);
        n.v = -51.0;
        n.eta = 0.3;
        n.theta = -90.0;
        n.lambda_0 = 1.0e9;
        assert_eq!(n.step(0.0), 1);
        assert!((n.v - n.v_reset).abs() < 1e-12);
        assert!((n.eta - 5.298503743757805).abs() < 1e-15);
    }
    #[test]
    fn gif_pop_invalid_input_preserves_state() {
        let mut n = GIFPopulationNeuron::new(42);
        n.v = -62.0;
        n.eta = 0.75;
        let before = (n.v, n.eta);
        assert_eq!(n.step(f64::NAN), 0);
        n.tau_m = 0.0;
        assert_eq!(n.step(1.0), 0);
        assert_eq!((n.v, n.eta), before);
    }
    #[test]
    fn gif_pop_seeded_reset_replays() {
        let mut n = GIFPopulationNeuron::new(123);
        n.theta = -90.0;
        n.lambda_0 = 1.0e9;
        let first: Vec<i32> = (0..3).map(|_| n.step(0.0)).collect();
        n.reset();
        let replay: Vec<i32> = (0..3).map(|_| n.step(0.0)).collect();
        assert_eq!(first, replay);
        assert!(n.eta > n.eta_increment);
    }
    #[test]
    fn gif_pop_negative_drive_remains_finite() {
        let mut n = GIFPopulationNeuron::new(42);
        for _ in 0..200 {
            n.step(-30.0);
        }
        assert!(n.v.is_finite());
        assert!(n.eta.is_finite());
    }

    // -- AvRonCardiac --
    #[test]
    fn avron_rk4_reference_point() {
        let mut n = AvRonCardiacNeuron::new();
        n.v = -55.0;
        n.h = 0.55;
        n.n = 0.35;
        n.s = 0.45;
        assert_eq!(n.step(2.0), 0);
        assert!((n.v - (-50.0840498399381)).abs() < 1e-12);
        assert!((n.h - 0.5506609782132562).abs() < 1e-15);
        assert!((n.n - 0.34988677751350306).abs() < 1e-15);
        assert!((n.s - 0.4500091998827305).abs() < 1e-15);
    }
    #[test]
    fn avron_invalid_input_preserves_state() {
        let mut n = AvRonCardiacNeuron::new();
        n.v = -52.0;
        n.h = 0.4;
        n.n = 0.2;
        n.s = 0.8;
        let before = (n.v, n.h, n.n, n.s);
        assert_eq!(n.step(f64::NAN), 0);
        n.dt = 0.0;
        assert_eq!(n.step(1.0), 0);
        assert_eq!((n.v, n.h, n.n, n.s), before);
    }
    #[test]
    fn avron_rejects_corrupted_gate_state() {
        let mut n = AvRonCardiacNeuron::new();
        n.h = 1.2;
        let before = (n.v, n.h, n.n, n.s);
        assert_eq!(n.step(1.0), 0);
        assert_eq!((n.v, n.h, n.n, n.s), before);
    }
    #[test]
    fn avron_slow_current_remains_physical() {
        let mut n = AvRonCardiacNeuron::new();
        for _ in 0..1000 {
            n.step(5.0);
            assert!((0.0..=1.0).contains(&n.h));
            assert!((0.0..=1.0).contains(&n.n));
            assert!((0.0..=1.0).contains(&n.s));
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn avron_reset_clears_state() {
        let mut n = AvRonCardiacNeuron::new();
        n.v = -40.0;
        n.h = 0.1;
        n.n = 0.9;
        n.s = 0.2;
        n.reset();
        assert_eq!((n.v, n.h, n.n, n.s), (-60.0, 0.6, 0.3, 0.5));
    }

    // -- DurstewitzDopamine --
    #[test]
    fn durstewitz_low_activity_zero_input() {
        let mut n = DurstewitzDopamineNeuron::new();
        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        // NMDA tonic conductance can produce spontaneous activity
        assert!(n.v.is_finite());
    }
    #[test]
    fn durstewitz_reset_clears_state() {
        let mut n = DurstewitzDopamineNeuron::new();
        for _ in 0..100 {
            n.step(3.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn durstewitz_extreme_bounded() {
        let mut n = DurstewitzDopamineNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn durstewitz_d1_modulation() {
        // D1 dopamine should increase NMDA and shift Na activation
        let mut n_d1 = DurstewitzDopamineNeuron::new();
        n_d1.d1_level = 1.0;
        let mut n_no = DurstewitzDopamineNeuron::new();
        n_no.d1_level = 0.0;
        for _ in 0..1000 {
            n_d1.step(3.0);
        }
        for _ in 0..1000 {
            n_no.step(3.0);
        }
        // Both should remain stable; D1 changes effective conductances
        assert!(n_d1.v.is_finite() && n_no.v.is_finite());
    }
    #[test]
    fn durstewitz_mg_block() {
        let n = DurstewitzDopamineNeuron::new();
        // At rest (-65 mV), Mg²⁺ block should be high
        let block = 1.0 / (1.0 + n.mg * (-0.062 * n.v).exp() / 3.57);
        assert!(block < 0.1, "Mg²⁺ block at rest should be high: {}", block);
    }
    #[test]
    fn durstewitz_negative_no_crash() {
        let mut n = DurstewitzDopamineNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }

    // -- HillTononi --
    #[test]
    fn hill_tononi_silent_without_input() {
        let mut n = HillTononiNeuron::new();
        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        // May have some spontaneous activity due to Ih
        assert!(n.v.is_finite());
    }
    #[test]
    fn hill_tononi_reset_clears_state() {
        let mut n = HillTononiNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-65.0)).abs() < 1e-10);
        assert!((n.na_i - 5.0).abs() < 1e-10);
    }
    #[test]
    fn hill_tononi_extreme_bounded() {
        let mut n = HillTononiNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn hill_tononi_na_accumulation() {
        let mut n = HillTononiNeuron::new();
        for _ in 0..500 {
            n.step(5.0);
        }
        // Intracellular Na+ should remain finite and non-negative
        assert!(n.na_i.is_finite());
        assert!(n.na_i >= 0.0, "Na_i must be non-negative");
    }
    #[test]
    fn hill_tononi_negative_no_crash() {
        let mut n = HillTononiNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn hill_tononi_nan_no_panic() {
        let mut n = HillTononiNeuron::new();
        n.step(f64::NAN);
    }

    // -- BertramPhantom --
    #[test]
    fn bertram_silent_without_input() {
        let mut n = BertramPhantomBurster::new();
        for _ in 0..500 {
            n.step(0.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn bertram_reset_clears_state() {
        let mut n = BertramPhantomBurster::new();
        for _ in 0..1000 {
            n.step(200.0);
        }
        n.reset();
        assert!((n.v - (-50.0)).abs() < 1e-10);
        assert!((n.s1 - 0.1).abs() < 1e-10);
        assert!((n.s2 - 0.1).abs() < 1e-10);
    }
    #[test]
    fn bertram_extreme_bounded() {
        let mut n = BertramPhantomBurster::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn bertram_dual_slow_vars() {
        let mut n = BertramPhantomBurster::new();
        for _ in 0..10000 {
            n.step(200.0);
        }
        // Both slow variables should evolve
        assert!(n.s1 >= 0.0 && n.s1 <= 1.0, "s1={}", n.s1);
        assert!(n.s2 >= 0.0 && n.s2 <= 1.0, "s2={}", n.s2);
    }
    #[test]
    fn bertram_negative_no_crash() {
        let mut n = BertramPhantomBurster::new();
        for _ in 0..200 {
            n.step(-100.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn bertram_nan_no_panic() {
        let mut n = BertramPhantomBurster::new();
        n.step(f64::NAN);
    }

    // -- Yamada --
    #[test]
    fn yamada_silent_without_input() {
        let mut n = YamadaNeuron::new();
        let t: i32 = (0..500).map(|_| n.step(0.0)).sum();
        assert_eq!(t, 0);
    }
    #[test]
    fn yamada_reset_clears_state() {
        let mut n = YamadaNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-60.0)).abs() < 1e-10);
        assert!((n.q - 0.0).abs() < 1e-10);
    }
    #[test]
    fn yamada_extreme_bounded() {
        let mut n = YamadaNeuron::new();
        for _ in 0..200 {
            n.step(1e4);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn yamada_slow_q_adapts() {
        let mut n = YamadaNeuron::new();
        for _ in 0..2000 {
            n.step(5.0);
        }
        assert!(n.q > 0.0, "slow variable q should activate during spiking");
    }
    #[test]
    fn yamada_negative_no_crash() {
        let mut n = YamadaNeuron::new();
        for _ in 0..200 {
            n.step(-10.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn yamada_nan_no_panic() {
        let mut n = YamadaNeuron::new();
        n.step(f64::NAN);
    }
}