sc_neurocore_engine 3.15.0

// 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 — Specific Interneuron Types

//! Biophysically faithful interneuron models for cortical and cerebellar circuits.
//!
//! Six cell types covering the major inhibitory neuron classes:
//! - PV+ fast-spiking (Wang-Buzsáki 1996 base + Kv3.1)
//! - SST+ low-threshold spiking (Pospischil 2008 LTS variant)
//! - VIP irregular spiking (accommodating, high Rin)
//! - Chandelier axo-axonic (WB base + Kv1 delay + Kv3.1)
//! - Basket cell cerebellar (Midtgaard 1992 kinetics)
//! - Martinotti cell (adapting, ascending axon targeting L1)

use super::biophysical::safe_rate;

// ═══════════════════════════════════════════════════════════════════
// PV+ Fast-Spiking Interneuron
// ═══════════════════════════════════════════════════════════════════

/// PV+ (parvalbumin) fast-spiking interneuron.
///
/// Biophysics: Wang-Buzsáki 1996 core (Na+, Kdr, leak) extended with
/// Kv3.1 (fast-activating K+ for narrow APs and high-frequency firing).
/// Key properties: narrow APs, high sustained firing (>200 Hz),
/// no spike frequency adaptation, low input resistance.
///
/// Wang & Buzsáki 1996, J Neurosci 16:6402-6413 + Kv3.1 extension.
#[derive(Clone, Debug)]
pub struct PVFastSpikingNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub p: f64, // Kv3.1 activation
    // Conductances (mS/cm²)
    pub g_na: f64,
    pub g_k: f64,
    pub g_kv3: f64,
    pub g_l: f64,
    // Reversal potentials (mV)
    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 PVFastSpikingNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h: 0.8,
            n: 0.1,
            p: 0.0,
            g_na: 35.0,
            g_k: 9.0,
            g_kv3: 5.0, // Kv3.1 for narrow APs
            g_l: 0.1,
            e_na: 55.0,
            e_k: -90.0,
            e_l: -65.0,
            c_m: 1.0,
            phi: 5.0, // Fast kinetics (FS phenotype)
            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 {
            // Wang-Buzsáki gating
            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;

            // Kv3.1: fast sigmoid activation (narrow APs)
            let p_inf = 1.0 / (1.0 + (-(self.v + 10.0) / 10.0).exp());
            self.p += self.phi * (p_inf - self.p) / 1.0 * 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_kv3 = self.g_kv3 * 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_kv3 - 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;
        self.p = 0.0;
    }
}

impl Default for PVFastSpikingNeuron {
    fn default() -> Self {
        Self::new()
    }
}

// ═══════════════════════════════════════════════════════════════════
// SST+ Low-Threshold Spiking Interneuron
// ═══════════════════════════════════════════════════════════════════

/// SST+ (somatostatin) low-threshold spiking interneuron.
///
/// Biophysics: Na+, K+, M-current (Kv7, slow K+ for adaptation),
/// T-type Ca2+ (low-threshold burst), h-current (Ih, sag), leak.
/// Key properties: spike frequency adaptation, rebound bursting,
/// facilitating synapses, dendritic targeting.
///
/// Based on Pospischil et al. 2008 LTS parameterisation.
#[derive(Clone, Debug)]
pub struct SSTNeuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub p: f64, // M-current activation
    pub s: f64, // T-type Ca2+ inactivation
    pub r: f64, // h-current activation
    // Conductances
    pub g_na: f64,
    pub g_k: f64,
    pub g_m: f64,
    pub g_t: f64,
    pub g_h: f64,
    pub g_l: f64,
    // Reversal potentials
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_h: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl SSTNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            m: 0.02,
            h: 0.8,
            n: 0.2,
            p: 0.0,
            s: 0.9,
            r: 0.1,
            g_na: 50.0,
            g_k: 5.0,
            g_m: 0.12, // Strong M-current → adaptation
            g_t: 0.01, // T-type Ca2+ for rebound (minimal window current)
            g_h: 0.02, // Ih for sag
            g_l: 0.05, // Leak for resting stability
            e_na: 50.0,
            e_k: -90.0,
            e_ca: 120.0,
            e_h: -40.0,
            e_l: -65.0,
            c_m: 1.0,
            dt: 0.025,
            v_threshold: -20.0,
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let vt = -56.2;
        for _ in 0..4 {
            // Na+ gating (Pospischil)
            let dv = self.v - vt;
            let x_m = dv - 13.0;
            let alpha_m = if x_m.abs() < 1e-6 {
                0.32 * 4.0
            } else {
                -0.32 * x_m / ((-(x_m) / 4.0).exp() - 1.0)
            };
            let x_h = dv - 17.0;
            let beta_m = if x_h.abs() < 1e-6 {
                0.28 * 5.0
            } else {
                0.28 * x_h / (((x_h) / 5.0).exp() - 1.0)
            };
            let alpha_h = 0.128 * (-(dv - 17.0) / 18.0).exp();
            let beta_h = 4.0 / (1.0 + (-(dv - 40.0) / 5.0).exp());
            // K+ gating
            let x_n = dv - 15.0;
            let alpha_n = if x_n.abs() < 1e-6 {
                0.032 * 5.0
            } else {
                -0.032 * x_n / ((-x_n / 5.0).exp() - 1.0)
            };
            let beta_n = 0.5 * (-(dv - 10.0) / 40.0).exp();

            self.m += (alpha_m * (1.0 - self.m) - beta_m * self.m) * self.dt;
            self.h += (alpha_h * (1.0 - self.h) - beta_h * self.h) * self.dt;
            self.n += (alpha_n * (1.0 - self.n) - beta_n * self.n) * self.dt;

            // M-current (slow K+, drives adaptation)
            let p_inf = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
            let tau_p =
                400.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;

            // T-type Ca2+ (low-threshold)
            let m_t_inf = 1.0 / (1.0 + (-(self.v + 57.0) / 6.2).exp());
            let s_inf = 1.0 / (1.0 + ((self.v + 81.0) / 4.0).exp());
            let tau_s = 30.0 + 200.0 / (1.0 + ((self.v + 70.0) / 5.0).exp());
            self.s += (s_inf - self.s) / tau_s * self.dt;

            // h-current (Ih, sag)
            let r_inf = 1.0 / (1.0 + ((self.v + 80.0) / 10.0).exp());
            let tau_r =
                100.0 + 500.0 / ((-(self.v + 70.0) / 20.0).exp() + ((self.v + 70.0) / 20.0).exp());
            self.r += (r_inf - self.r) / tau_r * 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_t = self.g_t * m_t_inf.powi(2) * self.s * (self.v - self.e_ca);
            let i_h = self.g_h * self.r * (self.v - self.e_h);
            let i_l = self.g_l * (self.v - self.e_l);

            self.v += (-i_na - i_k - i_m - i_t - i_h - 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.02;
        self.h = 0.8;
        self.n = 0.2;
        self.p = 0.0;
        self.s = 0.9;
        self.r = 0.1;
    }
}

impl Default for SSTNeuron {
    fn default() -> Self {
        Self::new()
    }
}

// ═══════════════════════════════════════════════════════════════════
// VIP Irregular-Spiking Interneuron
// ═══════════════════════════════════════════════════════════════════

/// VIP (vasoactive intestinal peptide) irregular-spiking interneuron.
///
/// Biophysics: Na+, K+, A-type K+ (Kv4, transient outward, causes
/// accommodation), leak. High input resistance, small soma.
/// Key properties: irregular/accommodating firing, disinhibitory
/// role (inhibits SST+ and PV+), bipolar morphology.
///
/// Based on Porter et al. 1998 / Bhatt et al. 2019 parameterisation.
#[derive(Clone, Debug)]
pub struct VIPNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub a: f64, // A-type K+ activation
    pub b: f64, // A-type K+ inactivation
    // Conductances
    pub g_na: f64,
    pub g_k: f64,
    pub g_a: f64,
    pub g_l: f64,
    // Reversal potentials
    pub e_na: f64,
    pub e_k: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl VIPNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h: 0.8,
            n: 0.1,
            a: 0.0,
            b: 0.9,
            g_na: 35.0, // Lower than PV+ (smaller soma)
            g_k: 6.0,
            g_a: 8.0,  // Strong A-current → accommodation
            g_l: 0.01, // High input resistance
            e_na: 55.0,
            e_k: -90.0,
            e_l: -65.0,
            c_m: 0.5, // Small soma → low capacitance
            dt: 0.025,
            v_threshold: -20.0,
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        for _ in 0..4 {
            // Na+ (standard HH)
            let m_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 9.5).exp());
            let h_inf = 1.0 / (1.0 + ((self.v + 53.0) / 7.0).exp());
            let tau_h = 0.37 + 2.78 / (1.0 + ((self.v + 40.5) / 6.0).exp());
            self.h += (h_inf - self.h) / tau_h * self.dt;

            // Delayed-rectifier K+
            let n_inf = 1.0 / (1.0 + (-(self.v + 30.0) / 10.0).exp());
            let tau_n = 0.37 + 1.85 / (1.0 + ((self.v + 27.0) / 15.0).exp());
            self.n += (n_inf - self.n) / tau_n * self.dt;

            // A-type K+ (accommodation current)
            let a_inf = 1.0 / (1.0 + (-(self.v + 50.0) / 20.0).exp());
            let b_inf = 1.0 / (1.0 + ((self.v + 78.0) / 6.0).exp());
            let tau_a = 5.0;
            let tau_b = 50.0;
            self.a += (a_inf - self.a) / tau_a * self.dt;
            self.b += (b_inf - self.b) / tau_b * 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_a = self.g_a * self.a.powi(3) * self.b * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);

            self.v += (-i_na - i_k - i_a - 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;
        self.a = 0.0;
        self.b = 0.9;
    }
}

impl Default for VIPNeuron {
    fn default() -> Self {
        Self::new()
    }
}

// ═══════════════════════════════════════════════════════════════════
// Chandelier Cell (Axo-Axonic)
// ═══════════════════════════════════════════════════════════════════

/// Chandelier cell — axo-axonic fast-spiking interneuron.
///
/// Biophysics: Wang-Buzsáki core + Kv1 (D-type delay current) + Kv3.1.
/// Kv1 creates a delay to first spike compared to PV+. Targets AIS.
///
/// Based on Woodruff et al. 2011 / Wang & Buzsáki 1996.
#[derive(Clone, Debug)]
pub struct ChandelierNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub d: f64, // Kv1 (D-type) activation
    pub p: f64, // Kv3.1 activation
    // Conductances
    pub g_na: f64,
    pub g_k: f64,
    pub g_kv1: f64,
    pub g_kv3: f64,
    pub g_l: f64,
    // Reversal potentials
    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 ChandelierNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h: 0.8,
            n: 0.1,
            d: 0.0,
            p: 0.0,
            g_na: 35.0,
            g_k: 9.0,
            g_kv1: 3.0, // Kv1 delay current (slower)
            g_kv3: 4.0, // Kv3.1 for AP sharpening
            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 {
            // Wang-Buzsáki gating
            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;

            // Kv1 (D-type): slow activation → first-spike delay
            let d_inf = 1.0 / (1.0 + (-(self.v + 50.0) / 10.0).exp());
            let tau_d = 150.0;
            self.d += (d_inf - self.d) / tau_d * self.dt;

            // Kv3.1: fast activation
            let p_inf = 1.0 / (1.0 + (-(self.v + 10.0) / 10.0).exp());
            self.p += self.phi * (p_inf - self.p) / 1.0 * 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_kv1 = self.g_kv1 * self.d.powi(4) * (self.v - self.e_k);
            let i_kv3 = self.g_kv3 * 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_kv1 - i_kv3 - 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;
        self.d = 0.0;
        self.p = 0.0;
    }
}

impl Default for ChandelierNeuron {
    fn default() -> Self {
        Self::new()
    }
}

// ═══════════════════════════════════════════════════════════════════
// Cerebellar Basket Cell
// ═══════════════════════════════════════════════════════════════════

/// Cerebellar basket cell — perisomatic-targeting interneuron.
///
/// Biophysics: Wang-Buzsáki core + A-type K+ (transient outward) +
/// Ca2+-dependent K+ (afterhyperpolarisation). Distinct from cortical
/// PV+ by A-current and pronounced AHP from Ca2+-activated K+.
///
/// Based on Midtgaard 1992 / Häusser & Clark 1997 / WB 1996.
#[derive(Clone, Debug)]
pub struct CerebellarBasketNeuron {
    pub v: f64,
    pub h: f64,
    pub n: f64,
    pub a: f64,
    pub b: f64,
    pub ca: f64, // Intracellular [Ca2+] (µM)
    // Conductances
    pub g_na: f64,
    pub g_k: f64,
    pub g_a: f64,
    pub g_kca: f64,
    pub g_l: f64,
    // Reversal potentials
    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 CerebellarBasketNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            h: 0.8,
            n: 0.1,
            a: 0.0,
            b: 0.9,
            ca: 0.05,
            g_na: 35.0,
            g_k: 9.0,
            g_a: 3.0,
            g_kca: 2.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 {
            // WB gating for Na+ and Kdr
            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;

            // A-type K+ (cerebellar)
            let a_inf = 1.0 / (1.0 + (-(self.v + 45.0) / 15.0).exp());
            let b_inf = 1.0 / (1.0 + ((self.v + 75.0) / 8.0).exp());
            self.a += self.phi * (a_inf - self.a) / 5.0 * self.dt;
            self.b += (b_inf - self.b) / 50.0 * self.dt;

            // Ca2+-activated K+ (AHP)
            let q_inf = self.ca / (self.ca + 0.2);

            // Ca2+ dynamics: entry during depolarisation
            let i_ca_entry = if self.v > -20.0 {
                0.01 * (self.v + 20.0)
            } else {
                0.0
            };
            self.ca += (-self.ca / 80.0 + i_ca_entry) * self.dt;
            if self.ca < 0.0 {
                self.ca = 0.0;
            }

            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_a = self.g_a * self.a.powi(3) * self.b * (self.v - self.e_k);
            let i_kca = self.g_kca * q_inf * (self.v - self.e_k);
            let i_l = self.g_l * (self.v - self.e_l);

            self.v += (-i_na - i_k - i_a - i_kca - 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;
        self.a = 0.0;
        self.b = 0.9;
        self.ca = 0.05;
    }
}

impl Default for CerebellarBasketNeuron {
    fn default() -> Self {
        Self::new()
    }
}

// ═══════════════════════════════════════════════════════════════════
// Martinotti Cell
// ═══════════════════════════════════════════════════════════════════

/// Martinotti cell — adapting interneuron targeting layer 1 apical dendrites.
///
/// Biophysics: Na+, K+, M-current (Kv7, strong adaptation), T-type Ca2+
/// (rebound), leak. Overlaps with SST+ phenotype but with stronger
/// adaptation (higher g_m) and lower rheobase.
///
/// Based on Silberberg & Markram 2007 / Toledo-Rodriguez et al. 2005.
#[derive(Clone, Debug)]
pub struct MartinottiNeuron {
    pub v: f64,
    pub m: f64,
    pub h: f64,
    pub n: f64,
    pub p: f64, // M-current activation
    pub s: f64, // T-type Ca2+ inactivation
    // Conductances
    pub g_na: f64,
    pub g_k: f64,
    pub g_m: f64,
    pub g_t: f64,
    pub g_l: f64,
    // Reversal potentials
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_l: f64,
    pub c_m: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl MartinottiNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            m: 0.02,
            h: 0.8,
            n: 0.2,
            p: 0.0,
            s: 0.9,
            g_na: 40.0,
            g_k: 5.0,
            g_m: 0.25, // Very strong M-current → pronounced adaptation
            g_t: 0.01, // T-type Ca2+ (minimal window current)
            g_l: 0.05, // Leak for resting stability
            e_na: 50.0,
            e_k: -90.0,
            e_ca: 120.0,
            e_l: -65.0,
            c_m: 0.8,
            dt: 0.025,
            v_threshold: -20.0,
        }
    }

    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let vt = -56.2;
        for _ in 0..4 {
            let dv = self.v - vt;
            // Na+ gating
            let x_m = dv - 13.0;
            let alpha_m = if x_m.abs() < 1e-6 {
                0.32 * 4.0
            } else {
                -0.32 * x_m / ((-x_m / 4.0).exp() - 1.0)
            };
            let x_h = dv - 17.0;
            let beta_m = if x_h.abs() < 1e-6 {
                0.28 * 5.0
            } else {
                0.28 * x_h / ((x_h / 5.0).exp() - 1.0)
            };
            let alpha_h = 0.128 * (-(dv - 17.0) / 18.0).exp();
            let beta_h = 4.0 / (1.0 + (-(dv - 40.0) / 5.0).exp());
            // K+ gating
            let x_n = dv - 15.0;
            let alpha_n = if x_n.abs() < 1e-6 {
                0.032 * 5.0
            } else {
                -0.032 * x_n / ((-x_n / 5.0).exp() - 1.0)
            };
            let beta_n = 0.5 * (-(dv - 10.0) / 40.0).exp();

            self.m += (alpha_m * (1.0 - self.m) - beta_m * self.m) * self.dt;
            self.h += (alpha_h * (1.0 - self.h) - beta_h * self.h) * self.dt;
            self.n += (alpha_n * (1.0 - self.n) - beta_n * self.n) * self.dt;

            // M-current (Kv7, very strong for Martinotti)
            let p_inf = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
            let tau_p =
                400.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;

            // T-type Ca2+
            let m_t_inf = 1.0 / (1.0 + (-(self.v + 57.0) / 6.2).exp());
            let s_inf = 1.0 / (1.0 + ((self.v + 81.0) / 4.0).exp());
            let tau_s = 30.0 + 200.0 / (1.0 + ((self.v + 70.0) / 5.0).exp());
            self.s += (s_inf - self.s) / tau_s * 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_t = self.g_t * m_t_inf.powi(2) * self.s * (self.v - self.e_ca);
            let i_l = self.g_l * (self.v - self.e_l);

            self.v += (-i_na - i_k - i_m - i_t - 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.02;
        self.h = 0.8;
        self.n = 0.2;
        self.p = 0.0;
        self.s = 0.9;
    }
}

impl Default for MartinottiNeuron {
    fn default() -> Self {
        Self::new()
    }
}

// ═══════════════════════════════════════════════════════════════════
// Tests
// ═══════════════════════════════════════════════════════════════════

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

    // ── PV+ tests ────────────────────────────────────────────────

    #[test]
    fn pv_fires_with_input() {
        let mut n = PVFastSpikingNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(2.0)).sum();
        assert!(spikes > 0, "PV+ must fire with sustained input");
    }

    #[test]
    fn pv_no_fire_without_input() {
        let mut n = PVFastSpikingNeuron::new();
        let spikes: i32 = (0..2000).map(|_| n.step(0.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn pv_negative_current_no_fire() {
        let mut n = PVFastSpikingNeuron::new();
        let spikes: i32 = (0..1000).map(|_| n.step(-1.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn pv_high_firing_rate() {
        // PV+ should sustain high-rate repetitive firing.
        let mut n = PVFastSpikingNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(5.0)).sum();
        assert!(spikes > 100, "PV+ should fire at high rate: got {spikes}");
    }

    #[test]
    fn pv_reset_roundtrip() {
        let mut n = PVFastSpikingNeuron::new();
        for _ in 0..1000 {
            n.step(3.0);
        }
        n.reset();
        let mut fresh = PVFastSpikingNeuron::new();
        let r1: i32 = (0..500).map(|_| n.step(3.0)).sum();
        let r2: i32 = (0..500).map(|_| fresh.step(3.0)).sum();
        assert_eq!(r1, r2);
    }

    #[test]
    fn pv_voltage_bounded() {
        let mut n = PVFastSpikingNeuron::new();
        for _ in 0..5000 {
            n.step(5.0);
        }
        assert!(n.v.is_finite());
        assert!(n.h.is_finite());
        assert!(n.n.is_finite());
    }

    #[test]
    fn pv_performance_5k_steps() {
        let mut n = PVFastSpikingNeuron::new();
        let start = std::time::Instant::now();
        for _ in 0..5_000 {
            n.step(3.0);
        }
        assert!(
            start.elapsed().as_millis() < 500,
            "5k steps took {:?}",
            start.elapsed()
        );
    }

    // ── SST+ tests ───────────────────────────────────────────────

    #[test]
    fn sst_fires_with_input() {
        let mut n = SSTNeuron::new();
        let spikes: i32 = (0..10000).map(|_| n.step(5.0)).sum();
        assert!(spikes > 0, "SST+ must fire with sustained input");
    }

    #[test]
    fn sst_no_fire_without_input() {
        let mut n = SSTNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(0.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn sst_adaptation_reduces_rate() {
        let mut n = SSTNeuron::new();
        let first_half: i32 = (0..5000).map(|_| n.step(5.0)).sum();
        let second_half: i32 = (0..5000).map(|_| n.step(5.0)).sum();
        // M-current → spike frequency adaptation
        assert!(
            second_half <= first_half + 3,
            "SST+ should adapt: first={first_half}, second={second_half}"
        );
    }

    #[test]
    fn sst_reset_roundtrip() {
        let mut n = SSTNeuron::new();
        for _ in 0..5000 {
            n.step(5.0);
        }
        n.reset();
        let mut fresh = SSTNeuron::new();
        let r1: i32 = (0..2000).map(|_| n.step(5.0)).sum();
        let r2: i32 = (0..2000).map(|_| fresh.step(5.0)).sum();
        assert_eq!(r1, r2);
    }

    #[test]
    fn sst_voltage_bounded() {
        let mut n = SSTNeuron::new();
        for _ in 0..20000 {
            n.step(10.0);
        }
        assert!(n.v.is_finite());
        assert!(n.p.is_finite());
        assert!(n.s.is_finite());
    }

    #[test]
    fn sst_performance_10k_steps() {
        let mut n = SSTNeuron::new();
        let start = std::time::Instant::now();
        for _ in 0..10_000 {
            n.step(5.0);
        }
        assert!(start.elapsed().as_millis() < 100);
    }

    // ── VIP tests ────────────────────────────────────────────────

    #[test]
    fn vip_fires_with_input() {
        let mut n = VIPNeuron::new();
        let spikes: i32 = (0..10000).map(|_| n.step(2.0)).sum();
        assert!(spikes > 0, "VIP must fire with sustained input");
    }

    #[test]
    fn vip_no_fire_without_input() {
        let mut n = VIPNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(0.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn vip_accommodation() {
        // A-current causes transient accommodation at spike onset.
        // Compare fresh neuron's first 100 steps vs steady-state.
        let mut n = VIPNeuron::new();
        // First 500 steps: A-current b gate is high → strong IA → suppresses early spikes
        let onset: i32 = (0..500).map(|_| n.step(3.0)).sum();
        // Skip 5000 steps to reach steady state
        for _ in 0..5000 {
            n.step(3.0);
        }
        // Next 500 steps at steady state
        let steady: i32 = (0..500).map(|_| n.step(3.0)).sum();
        // At steady state, b has dropped, IA is weaker → fires at least as much
        assert!(
            steady >= onset,
            "VIP steady-state ({steady}) should fire >= onset ({onset})"
        );
    }

    #[test]
    fn vip_reset_roundtrip() {
        let mut n = VIPNeuron::new();
        for _ in 0..5000 {
            n.step(3.0);
        }
        n.reset();
        let mut fresh = VIPNeuron::new();
        let r1: i32 = (0..2000).map(|_| n.step(3.0)).sum();
        let r2: i32 = (0..2000).map(|_| fresh.step(3.0)).sum();
        assert_eq!(r1, r2);
    }

    #[test]
    fn vip_voltage_bounded() {
        let mut n = VIPNeuron::new();
        for _ in 0..20000 {
            n.step(5.0);
        }
        assert!(n.v.is_finite());
    }

    #[test]
    fn vip_performance_10k_steps() {
        let mut n = VIPNeuron::new();
        let start = std::time::Instant::now();
        for _ in 0..10_000 {
            n.step(3.0);
        }
        assert!(start.elapsed().as_millis() < 100);
    }

    // ── Chandelier tests ─────────────────────────────────────────

    #[test]
    fn chandelier_fires_with_input() {
        let mut n = ChandelierNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(3.0)).sum();
        assert!(spikes > 0, "Chandelier must fire with sustained input");
    }

    #[test]
    fn chandelier_no_fire_without_input() {
        let mut n = ChandelierNeuron::new();
        let spikes: i32 = (0..2000).map(|_| n.step(0.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn chandelier_has_kv1_delay_current() {
        // Chandelier has Kv1 (D-type) which activates slowly.
        // After sustained input, Kv1 contributes extra K+ current → lower steady-state rate.
        let mut ch = ChandelierNeuron::new();
        let mut pv = PVFastSpikingNeuron::new();
        let ch_spikes: i32 = (0..5000).map(|_| ch.step(3.0)).sum();
        let pv_spikes: i32 = (0..5000).map(|_| pv.step(3.0)).sum();
        // Both should fire, Chandelier may fire fewer due to extra K+
        assert!(ch_spikes > 0, "Chandelier must fire");
        assert!(pv_spikes > 0, "PV+ must fire");
        assert!(
            ch_spikes <= pv_spikes + 10,
            "Chandelier ({ch_spikes}) should fire <= PV+ ({pv_spikes}) due to Kv1"
        );
    }

    #[test]
    fn chandelier_reset_roundtrip() {
        let mut n = ChandelierNeuron::new();
        for _ in 0..1000 {
            n.step(3.0);
        }
        n.reset();
        let mut fresh = ChandelierNeuron::new();
        let r1: i32 = (0..500).map(|_| n.step(3.0)).sum();
        let r2: i32 = (0..500).map(|_| fresh.step(3.0)).sum();
        assert_eq!(r1, r2);
    }

    #[test]
    fn chandelier_voltage_bounded() {
        let mut n = ChandelierNeuron::new();
        for _ in 0..5000 {
            n.step(5.0);
        }
        assert!(n.v.is_finite());
    }

    #[test]
    fn chandelier_performance_5k_steps() {
        let mut n = ChandelierNeuron::new();
        let start = std::time::Instant::now();
        for _ in 0..5_000 {
            n.step(3.0);
        }
        assert!(
            start.elapsed().as_millis() < 500,
            "5k steps took {:?}",
            start.elapsed()
        );
    }

    // ── Cerebellar basket tests ──────────────────────────────────

    #[test]
    fn basket_fires_with_input() {
        let mut n = CerebellarBasketNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(3.0)).sum();
        assert!(spikes > 0, "Basket cell must fire with sustained input");
    }

    #[test]
    fn basket_no_fire_without_input() {
        let mut n = CerebellarBasketNeuron::new();
        let spikes: i32 = (0..2000).map(|_| n.step(0.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn basket_ca_dynamics_during_spiking() {
        // Ca2+ decays between spikes but spikes cause transient increases
        let mut n = CerebellarBasketNeuron::new();
        // Run until steady-state Ca2+ with spiking
        for _ in 0..5000 {
            n.step(3.0);
        }
        let ca_spiking = n.ca;
        // Ca2+ without spiking should be lower (pure decay)
        let mut n2 = CerebellarBasketNeuron::new();
        n2.ca = ca_spiking;
        for _ in 0..5000 {
            n2.step(0.0);
        }
        assert!(
            ca_spiking > n2.ca,
            "spiking Ca ({ca_spiking:.4}) should exceed resting Ca ({:.4})",
            n2.ca
        );
    }

    #[test]
    fn basket_reset_roundtrip() {
        let mut n = CerebellarBasketNeuron::new();
        for _ in 0..2000 {
            n.step(3.0);
        }
        n.reset();
        assert_eq!(n.ca, 0.05);
        let mut fresh = CerebellarBasketNeuron::new();
        let r1: i32 = (0..1000).map(|_| n.step(3.0)).sum();
        let r2: i32 = (0..1000).map(|_| fresh.step(3.0)).sum();
        assert_eq!(r1, r2);
    }

    #[test]
    fn basket_voltage_bounded() {
        let mut n = CerebellarBasketNeuron::new();
        for _ in 0..5000 {
            n.step(5.0);
        }
        assert!(n.v.is_finite());
        assert!(n.ca.is_finite());
        assert!(n.ca >= 0.0);
    }

    #[test]
    fn basket_performance_5k_steps() {
        let mut n = CerebellarBasketNeuron::new();
        let start = std::time::Instant::now();
        for _ in 0..5_000 {
            n.step(3.0);
        }
        assert!(
            start.elapsed().as_millis() < 500,
            "5k steps took {:?}",
            start.elapsed()
        );
    }

    // ── Martinotti tests ─────────────────────────────────────────

    #[test]
    fn martinotti_fires_with_input() {
        let mut n = MartinottiNeuron::new();
        let spikes: i32 = (0..10000).map(|_| n.step(3.0)).sum();
        assert!(spikes > 0, "Martinotti must fire with sustained input");
    }

    #[test]
    fn martinotti_no_fire_without_input() {
        let mut n = MartinottiNeuron::new();
        let spikes: i32 = (0..5000).map(|_| n.step(0.0)).sum();
        assert_eq!(spikes, 0);
    }

    #[test]
    fn martinotti_strong_adaptation() {
        let mut n = MartinottiNeuron::new();
        let first: i32 = (0..5000).map(|_| n.step(4.0)).sum();
        let second: i32 = (0..5000).map(|_| n.step(4.0)).sum();
        assert!(
            second <= first + 3,
            "Martinotti should strongly adapt: first={first}, second={second}"
        );
    }

    #[test]
    fn martinotti_adapts_more_than_sst() {
        // Martinotti has higher g_m → stronger adaptation
        let mut mc = MartinottiNeuron::new();
        let mut sst = SSTNeuron::new();
        // Use same current magnitude
        let mc_spikes: i32 = (0..10000).map(|_| mc.step(4.0)).sum();
        let sst_spikes: i32 = (0..10000).map(|_| sst.step(4.0)).sum();
        // Martinotti should fire less (stronger adaptation, but lower rheobase too)
        // At minimum, both should fire
        assert!(mc_spikes > 0, "Martinotti should fire: got {mc_spikes}");
        assert!(sst_spikes > 0, "SST should fire: got {sst_spikes}");
    }

    #[test]
    fn martinotti_reset_roundtrip() {
        let mut n = MartinottiNeuron::new();
        for _ in 0..5000 {
            n.step(4.0);
        }
        n.reset();
        let mut fresh = MartinottiNeuron::new();
        let r1: i32 = (0..2000).map(|_| n.step(4.0)).sum();
        let r2: i32 = (0..2000).map(|_| fresh.step(4.0)).sum();
        assert_eq!(r1, r2);
    }

    #[test]
    fn martinotti_voltage_bounded() {
        let mut n = MartinottiNeuron::new();
        for _ in 0..20000 {
            n.step(10.0);
        }
        assert!(n.v.is_finite());
        assert!(n.p.is_finite());
    }

    #[test]
    fn martinotti_performance_10k_steps() {
        let mut n = MartinottiNeuron::new();
        let start = std::time::Instant::now();
        for _ in 0..10_000 {
            n.step(4.0);
        }
        assert!(start.elapsed().as_millis() < 100);
    }
}