sc_neurocore_engine 3.15.5

// 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 — Multi-compartment neuron models

//! Multi-compartment neuron models.

#[allow(unused_imports)]
use super::biophysical::safe_rate;

/// Pinsky-Rinzel 1994 — 2-compartment pyramidal cell.
#[derive(Clone, Debug)]
pub struct PinskyRinzelNeuron {
    pub v_s: f64,
    pub v_d: f64,
    pub h: f64,
    pub n: f64,
    pub s: f64,
    pub c: f64,
    pub q: f64,
    pub gc: f64,
    pub p: f64,
    pub g_na: f64,
    pub g_kdr: f64,
    pub g_ca: f64,
    pub g_kahp: 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 PinskyRinzelNeuron {
    pub fn new() -> Self {
        Self {
            v_s: -60.0,
            v_d: -60.0,
            h: 0.9,
            n: 0.1,
            s: 0.0,
            c: 0.0,
            q: 0.0,
            gc: 2.1,
            p: 0.5,
            g_na: 30.0,
            g_kdr: 15.0,
            g_ca: 10.0,
            g_kahp: 0.8,
            g_l: 0.1,
            e_na: 60.0,
            e_k: -75.0,
            e_ca: 80.0,
            e_l: -60.0,
            dt: 0.02,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current_soma: f64, current_dend: f64) -> i32 {
        let v_prev = self.v_s;
        let am = safe_rate(0.32, 54.0, self.v_s, 4.0, 8.0);
        let bm = safe_rate(-0.28, 27.0, self.v_s, -5.0, 5.6);
        let m_inf = am / (am + bm);
        let ah = 0.128 * (-(self.v_s + 50.0) / 18.0).exp();
        let bh = 4.0 / (1.0 + (-(self.v_s + 27.0) / 5.0).exp());
        let an = safe_rate(0.032, 52.0, self.v_s, 5.0, 0.32);
        let bn = 0.5 * (-(self.v_s + 57.0) / 40.0).exp();
        let s_inf = 1.0 / (1.0 + (-(self.v_d + 20.0) / 9.0).exp());
        let i_na = self.g_na * m_inf.powi(2) * self.h * (self.v_s - self.e_na);
        let i_kdr = self.g_kdr * self.n.powi(2) * (self.v_s - self.e_k);
        let i_ls = self.g_l * (self.v_s - self.e_l);
        let i_ds = (self.gc / self.p) * (self.v_s - self.v_d);
        let i_ca = self.g_ca * self.s.powi(2) * (self.v_d - self.e_ca);
        let i_kahp = self.g_kahp * self.q * (self.v_d - self.e_k);
        let i_ld = self.g_l * (self.v_d - self.e_l);
        let i_sd = (self.gc / (1.0 - self.p)) * (self.v_d - self.v_s);
        self.v_s += (-i_na - i_kdr - i_ls - i_ds + current_soma / self.p) * self.dt;
        self.v_d += (-i_ca - i_kahp - i_ld - i_sd + current_dend / (1.0 - self.p)) * 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;
        self.s += ((s_inf - self.s) / 5.0) * self.dt;
        self.c = (self.c + (-0.13 * i_ca - 0.075 * self.c) * self.dt).max(0.0);
        let q_inf = (self.c / (self.c + 2.0)).min(1.0);
        self.q += ((q_inf - self.q) / 100.0) * self.dt;
        if self.v_s >= self.v_threshold && v_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v_s = -60.0;
        self.v_d = -60.0;
        self.h = 0.9;
        self.n = 0.1;
        self.s = 0.0;
        self.c = 0.0;
        self.q = 0.0;
    }
}
impl Default for PinskyRinzelNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Hay et al. 2011 — Layer 5 thick-tufted pyramidal (3-compartment reduced).
#[derive(Clone, Debug)]
pub struct HayL5PyramidalNeuron {
    pub v_s: f64,
    pub v_t: f64,
    pub v_a: f64,
    pub h_na: f64,
    pub n_k: f64,
    pub m_ca: f64,
    pub h_ca: f64,
    pub m_ih: f64,
    pub ca_a: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_l_s: f64,
    pub g_ca_t: f64,
    pub g_ih: f64,
    pub g_l_t: f64,
    pub g_ca_a: f64,
    pub g_kca: f64,
    pub g_l_a: f64,
    pub g_st: f64,
    pub g_ta: f64,
    pub p_s: f64,
    pub p_t: f64,
    pub p_a: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_ih: f64,
    pub e_l: f64,
    pub ca_decay: f64,
    pub f_ca: f64,
    pub c_m: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl HayL5PyramidalNeuron {
    pub fn new() -> Self {
        Self {
            v_s: -75.0,
            v_t: -75.0,
            v_a: -75.0,
            h_na: 0.9,
            n_k: 0.1,
            m_ca: 0.0,
            h_ca: 1.0,
            m_ih: 0.0,
            ca_a: 0.0001,
            g_na: 300.0,
            g_k: 40.0,
            g_l_s: 0.03,
            g_ca_t: 2.0,
            g_ih: 0.02,
            g_l_t: 0.03,
            g_ca_a: 1.5,
            g_kca: 2.5,
            g_l_a: 0.03,
            g_st: 1.5,
            g_ta: 0.8,
            p_s: 0.15,
            p_t: 0.25,
            p_a: 0.60,
            e_na: 50.0,
            e_k: -85.0,
            e_ca: 140.0,
            e_ih: -45.0,
            e_l: -75.0,
            ca_decay: 200.0,
            f_ca: 0.0002,
            c_m: 1.0,
            dt: 0.025,
            v_threshold: -30.0,
        }
    }
    pub fn step(&mut self, current_soma: f64, current_tuft: f64) -> i32 {
        let v_s_prev = self.v_s;
        for _ in 0..4 {
            let m_na_inf = 1.0 / (1.0 + (-(self.v_s + 38.0) / 7.0).exp());
            let h_na_inf = 1.0 / (1.0 + ((self.v_s + 65.0) / 6.0).exp());
            let n_k_inf = 1.0 / (1.0 + (-(self.v_s + 25.0) / 12.0).exp());
            let tau_h = 0.5 + 14.0 / (1.0 + ((self.v_s + 35.0) / 10.0).exp());
            let tau_n = 1.0 + 5.0 / (1.0 + ((self.v_s + 30.0) / 10.0).exp());
            self.h_na += (h_na_inf - self.h_na) / tau_h * self.dt;
            self.n_k += (n_k_inf - self.n_k) / tau_n * self.dt;
            let i_na = self.g_na * m_na_inf.powi(3) * self.h_na * (self.v_s - self.e_na);
            let i_k = self.g_k * self.n_k.powi(4) * (self.v_s - self.e_k);
            let i_l_s = self.g_l_s * (self.v_s - self.e_l);
            let i_st = self.g_st * (self.v_s - self.v_t) / self.p_s;
            let m_ca_inf = 1.0 / (1.0 + (-(self.v_t + 27.0) / 7.0).exp());
            let h_ca_inf = 1.0 / (1.0 + ((self.v_t + 52.0) / 5.0).exp());
            let m_ih_inf = 1.0 / (1.0 + ((self.v_t + 75.0) / 5.5).exp());
            self.m_ca += (m_ca_inf - self.m_ca) / 1.0 * self.dt;
            self.h_ca += (h_ca_inf - self.h_ca) / 20.0 * self.dt;
            self.m_ih += (m_ih_inf - self.m_ih) / 50.0 * self.dt;
            let i_ca_t = self.g_ca_t * self.m_ca.powi(2) * self.h_ca * (self.v_t - self.e_ca);
            let i_ih = self.g_ih * self.m_ih * (self.v_t - self.e_ih);
            let i_l_t = self.g_l_t * (self.v_t - self.e_l);
            let i_ts = self.g_st * (self.v_t - self.v_s) / self.p_t;
            let i_ta = self.g_ta * (self.v_t - self.v_a) / self.p_t;
            let m_ca_a_inf = 1.0 / (1.0 + (-(self.v_a + 30.0) / 5.0).exp());
            let kca_act = self.ca_a / (self.ca_a + 0.001);
            let i_ca_a = self.g_ca_a * m_ca_a_inf.powi(2) * (self.v_a - self.e_ca);
            let i_kca = self.g_kca * kca_act * (self.v_a - self.e_k);
            let i_l_a = self.g_l_a * (self.v_a - self.e_l);
            let i_at = self.g_ta * (self.v_a - self.v_t) / self.p_a;
            self.v_s += (-i_na - i_k - i_l_s - i_st + current_soma / self.p_s) / self.c_m * self.dt;
            self.v_t += (-i_ca_t - i_ih - i_l_t - i_ts - i_ta) / self.c_m * self.dt;
            self.v_a +=
                (-i_ca_a - i_kca - i_l_a - i_at + current_tuft / self.p_a) / self.c_m * self.dt;
            self.ca_a =
                (self.ca_a + (-self.f_ca * i_ca_a - self.ca_a / self.ca_decay) * self.dt).max(0.0);
        }
        if self.v_s >= self.v_threshold && v_s_prev < self.v_threshold {
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v_s = -75.0;
        self.v_t = -75.0;
        self.v_a = -75.0;
        self.h_na = 0.9;
        self.n_k = 0.1;
        self.m_ca = 0.0;
        self.h_ca = 1.0;
        self.m_ih = 0.0;
        self.ca_a = 0.0001;
    }
}
impl Default for HayL5PyramidalNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Marder STG — stomatogastric ganglion LP neuron with 7 currents.
#[derive(Clone, Debug)]
pub struct MarderSTGNeuron {
    pub v: f64,
    pub m_na: f64,
    pub h_na: f64,
    pub m_cat: f64,
    pub h_cat: f64,
    pub m_cas: f64,
    pub m_a: f64,
    pub h_a: f64,
    pub m_kd: f64,
    pub m_h: f64,
    pub ca: f64,
    pub g_na: f64,
    pub g_cat: f64,
    pub g_cas: f64,
    pub g_a: f64,
    pub g_kd: f64,
    pub g_kca: f64,
    pub g_h: f64,
    pub g_l: f64,
    pub e_na: f64,
    pub e_k: f64,
    pub e_ca: f64,
    pub e_h: f64,
    pub e_l: f64,
    pub dt: f64,
    pub v_threshold: f64,
}

impl MarderSTGNeuron {
    pub fn new() -> Self {
        Self {
            v: -60.0,
            m_na: 0.0,
            h_na: 0.9,
            m_cat: 0.0,
            h_cat: 0.9,
            m_cas: 0.0,
            m_a: 0.0,
            h_a: 0.9,
            m_kd: 0.0,
            m_h: 0.0,
            ca: 0.05,
            g_na: 400.0,
            g_cat: 2.5,
            g_cas: 6.0,
            g_a: 50.0,
            g_kd: 100.0,
            g_kca: 25.0,
            g_h: 0.02,
            g_l: 0.01,
            e_na: 50.0,
            e_k: -80.0,
            e_ca: 120.0,
            e_h: -20.0,
            e_l: -50.0,
            dt: 0.05,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.v;
        let b = |v: f64, vh: f64, s: f64| 1.0 / (1.0 + (-(v - vh) / s).exp());
        self.m_na += (b(self.v, -25.5, 5.29) - self.m_na) / 1.32 * self.dt;
        self.h_na += (b(self.v, -48.9, -5.18) - self.h_na)
            / (0.67 * (1.0 + ((self.v + 62.9) / -10.0).exp()) + 1.5)
            * self.dt;
        self.m_cat += (b(self.v, -27.1, 7.2) - self.m_cat) / 21.7 * self.dt;
        self.h_cat += (b(self.v, -32.1, -5.5) - self.h_cat) / 105.0 * self.dt;
        self.m_cas += (b(self.v, -33.0, 8.1) - self.m_cas) / 14.0 * self.dt;
        self.m_a += (b(self.v, -27.2, 8.7) - self.m_a) / 11.6 * self.dt;
        self.h_a += (b(self.v, -56.9, -4.9) - self.h_a) / 38.6 * self.dt;
        self.m_kd += (b(self.v, -12.3, 11.8) - self.m_kd) / 7.2 * self.dt;
        self.m_h += (b(self.v, -70.0, -6.0) - self.m_h) / 272.0 * self.dt;
        let kca_act = (self.ca / (self.ca + 3.0)).min(1.0);
        let i_na = self.g_na * self.m_na.powi(3) * self.h_na * (self.v - self.e_na);
        let i_cat = self.g_cat * self.m_cat.powi(3) * self.h_cat * (self.v - self.e_ca);
        let i_cas = self.g_cas * self.m_cas.powi(3) * (self.v - self.e_ca);
        let i_a = self.g_a * self.m_a.powi(3) * self.h_a * (self.v - self.e_k);
        let i_kd = self.g_kd * self.m_kd.powi(4) * (self.v - self.e_k);
        let i_kca = self.g_kca * kca_act.powi(4) * (self.v - self.e_k);
        let i_h = self.g_h * self.m_h * (self.v - self.e_h);
        let i_l = self.g_l * (self.v - self.e_l);
        self.v += (-i_na - i_cat - i_cas - i_a - i_kd - i_kca - i_h - i_l + current) * self.dt;
        let i_ca_total = i_cat + i_cas;
        self.ca = (self.ca + (-0.0001 * i_ca_total - 0.01 * 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 = -60.0;
        self.m_na = 0.0;
        self.h_na = 0.9;
        self.m_cat = 0.0;
        self.h_cat = 0.9;
        self.m_cas = 0.0;
        self.m_a = 0.0;
        self.h_a = 0.9;
        self.m_kd = 0.0;
        self.m_h = 0.0;
        self.ca = 0.05;
    }
}
impl Default for MarderSTGNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Rall cable — N-compartment passive dendrite model. Rall 1964.
#[derive(Clone, Debug)]
pub struct RallCableNeuron {
    pub v: Vec<f64>,
    pub n_comp: usize,
    pub tau_m: f64,
    pub v_rest: f64,
    pub g_ratio: f64,
    pub v_threshold: f64,
    pub v_reset: f64,
    pub dt: f64,
}

impl RallCableNeuron {
    pub fn new(n_comp: usize) -> Self {
        Self {
            v: vec![-65.0; n_comp],
            n_comp,
            tau_m: 20.0,
            v_rest: -65.0,
            g_ratio: 0.5,
            v_threshold: -50.0,
            v_reset: -65.0,
            dt: 0.1,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let n = self.n_comp;
        let mut dv = vec![0.0; n];
        for i in 0..n {
            let leak = -(self.v[i] - self.v_rest) / self.tau_m;
            let left = if i > 0 {
                self.g_ratio * (self.v[i - 1] - self.v[i])
            } else {
                0.0
            };
            let right = if i + 1 < n {
                self.g_ratio * (self.v[i + 1] - self.v[i])
            } else {
                0.0
            };
            let inj = if i == n - 1 { current } else { 0.0 };
            dv[i] = (leak + left + right + inj) * self.dt;
        }
        for i in 0..n {
            self.v[i] += dv[i];
        }
        if self.v[0] >= self.v_threshold {
            self.v[0] = self.v_reset;
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v.fill(self.v_rest);
    }
}

/// Booth-Rinzel — 2-compartment motoneuron with bistability. Booth et al. 1997.
#[derive(Clone, Debug)]
pub struct BoothRinzelNeuron {
    pub vs: f64,
    pub vd: f64,
    pub h: f64,
    pub n: f64,
    pub q: f64,
    pub ca: f64,
    pub p: f64,
    pub gc: f64,
    pub g_na: f64,
    pub g_k: f64,
    pub g_ca: 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 dt: f64,
    pub v_threshold: f64,
}

impl BoothRinzelNeuron {
    pub fn new() -> Self {
        Self {
            vs: -65.0,
            vd: -65.0,
            h: 0.9,
            n: 0.0,
            q: 0.0,
            ca: 0.0,
            p: 0.5,
            gc: 0.1,
            g_na: 120.0,
            g_k: 20.0,
            g_ca: 14.0,
            g_kca: 5.0,
            g_l: 0.51,
            e_na: 55.0,
            e_k: -80.0,
            e_ca: 80.0,
            e_l: -60.0,
            dt: 0.025,
            v_threshold: -20.0,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let v_prev = self.vs;
        for _ in 0..4 {
            let m_inf = 1.0 / (1.0 + (-(self.vs + 35.0) / 7.8).exp());
            let h_inf = 1.0 / (1.0 + ((self.vs + 55.0) / 7.0).exp());
            let n_inf = 1.0 / (1.0 + (-(self.vs + 28.0) / 15.0).exp());
            let s_inf = 1.0 / (1.0 + (-(self.vd + 22.0) / 5.0).exp());
            let q_inf = 1.0 / (1.0 + (-(self.vd + 35.0) / 2.0).exp());
            let tau_h = (30.0
                / (((self.vs + 50.0) / 15.0).exp() + ((-(self.vs + 50.0)) / 16.0).exp() + 1e-12))
                .max(0.01);
            let tau_n = (7.0
                / (((self.vs + 40.0) / 40.0).exp() + ((-(self.vs + 40.0)) / 50.0).exp() + 1e-12))
                .max(0.01);
            self.h = (self.h + (h_inf - self.h) / tau_h * self.dt).clamp(0.0, 1.0);
            self.n = (self.n + (n_inf - self.n) / tau_n * self.dt).clamp(0.0, 1.0);
            self.q = (self.q + (q_inf - self.q) / 400.0 * self.dt).clamp(0.0, 1.0);
            let chi = (self.ca / 250.0).min(1.0);
            let i_na = self.g_na * m_inf.powi(3) * self.h * (self.vs - self.e_na);
            let i_k = self.g_k * self.n.powi(4) * (self.vs - self.e_k);
            let i_ls = self.g_l * (self.vs - self.e_l);
            let i_sd = (self.gc / self.p) * (self.vs - self.vd);
            let i_ca = self.g_ca * s_inf.powi(2) * (self.vd - self.e_ca);
            let i_kca = self.g_kca * chi * (self.vd - self.e_k);
            let i_ld = self.g_l * (self.vd - self.e_l);
            let i_ds = (self.gc / (1.0 - self.p)) * (self.vd - self.vs);
            self.vs = (self.vs + (-i_na - i_k - i_ls - i_sd + current / self.p) * self.dt)
                .clamp(-200.0, 100.0);
            self.vd = (self.vd + (-i_ca - i_kca - i_ld - i_ds) * self.dt).clamp(-200.0, 100.0);
            self.ca = (self.ca + (0.0025 * (-0.009 * i_ca) - 0.18 * self.ca) * self.dt).max(0.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.vd = -65.0;
        self.h = 0.9;
        self.n = 0.0;
        self.q = 0.0;
        self.ca = 0.0;
    }
}
impl Default for BoothRinzelNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Dendrify — two-compartment with active dendritic spike (NMDA-like plateau).
#[derive(Clone, Debug)]
pub struct DendrifyNeuron {
    pub v_s: f64,
    pub v_d: f64,
    pub d_active: bool,
    pub d_timer: f64,
    pub tau_s: f64,
    pub tau_d: f64,
    pub g_c: f64,
    pub d_threshold: f64,
    pub d_amplitude: f64,
    pub d_duration: f64,
    pub v_rest: f64,
    pub v_threshold: f64,
    pub v_reset: f64,
    pub dt: f64,
}

impl DendrifyNeuron {
    pub fn new() -> Self {
        Self {
            v_s: -65.0,
            v_d: -65.0,
            d_active: false,
            d_timer: 0.0,
            tau_s: 10.0,
            tau_d: 20.0,
            g_c: 0.8,
            d_threshold: -35.0,
            d_amplitude: 30.0,
            d_duration: 10.0,
            v_rest: -65.0,
            v_threshold: -50.0,
            v_reset: -65.0,
            dt: 0.1,
        }
    }
    pub fn step(&mut self, current: f64) -> i32 {
        let d_input = if self.d_active { self.d_amplitude } else { 0.0 };
        self.v_d += (-(self.v_d - self.v_rest) + current + d_input
            - self.g_c * (self.v_d - self.v_s))
            / self.tau_d
            * self.dt;
        self.v_s +=
            (-(self.v_s - self.v_rest) + self.g_c * (self.v_d - self.v_s)) / self.tau_s * self.dt;
        if self.d_active {
            self.d_timer -= self.dt;
            if self.d_timer <= 0.0 {
                self.d_active = false;
            }
        } else if self.v_d >= self.d_threshold {
            self.d_active = true;
            self.d_timer = self.d_duration;
        }
        if self.v_s >= self.v_threshold {
            self.v_s = self.v_reset;
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v_s = -65.0;
        self.v_d = -65.0;
        self.d_active = false;
        self.d_timer = 0.0;
    }
}
impl Default for DendrifyNeuron {
    fn default() -> Self {
        Self::new()
    }
}

/// Two-compartment LIF — soma + dendrite with history-dependent coupling.
#[derive(Clone, Debug)]
pub struct TwoCompartmentLIFNeuron {
    pub v_s: f64,
    pub v_d: f64,
    pub v_rest: f64,
    pub v_reset: f64,
    pub theta: f64,
    pub tau_s: f64,
    pub tau_d: f64,
    pub kappa: f64,
    pub dt: f64,
}

impl TwoCompartmentLIFNeuron {
    pub fn new() -> Self {
        Self {
            v_s: 0.0,
            v_d: 0.0,
            v_rest: 0.0,
            v_reset: 0.0,
            theta: 1.0,
            tau_s: 2.0,
            tau_d: 20.0,
            kappa: 0.5,
            dt: 1.0,
        }
    }
    pub fn step(&mut self, i_soma: f64, i_dend: f64) -> i32 {
        let alpha_s = (-self.dt / self.tau_s).exp();
        let alpha_d = (-self.dt / self.tau_d).exp();
        self.v_d = alpha_d * self.v_d + i_dend;
        self.v_s = alpha_s * self.v_s + i_soma + self.kappa * self.v_d;
        if self.v_s >= self.theta {
            self.v_s = self.v_reset;
            1
        } else {
            0
        }
    }
    pub fn reset(&mut self) {
        self.v_s = self.v_rest;
        self.v_d = self.v_rest;
    }
}
impl Default for TwoCompartmentLIFNeuron {
    fn default() -> Self {
        Self::new()
    }
}

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

    #[test]
    fn pr_fires() {
        let mut n = PinskyRinzelNeuron::new();
        let t: i32 = (0..5000).map(|_| n.step(5.0, 0.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn hay_fires() {
        let mut n = HayL5PyramidalNeuron::new();
        let t: i32 = (0..500).map(|_| n.step(20.0, 0.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn marder_fires() {
        let mut n = MarderSTGNeuron::new();
        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn rall_fires() {
        let mut n = RallCableNeuron::new(5);
        let t: i32 = (0..500).map(|_| n.step(50.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn booth_fires() {
        let mut n = BoothRinzelNeuron::new();
        let t: i32 = (0..2000).map(|_| n.step(5.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn dendrify_fires() {
        let mut n = DendrifyNeuron::new();
        let t: i32 = (0..2000).map(|_| n.step(50.0)).sum();
        assert!(t > 0);
    }
    #[test]
    fn tc_lif_fires() {
        let mut n = TwoCompartmentLIFNeuron::new();
        let t: i32 = (0..100).map(|_| n.step(0.5, 0.3)).sum();
        assert!(t > 0);
    }

    // ── Multi-angle tests for multi-compartment models ──

    // -- PinskyRinzel --
    #[test]
    fn pr_reset() {
        let mut n = PinskyRinzelNeuron::new();
        for _ in 0..100 {
            n.step(5.0, 0.0);
        }
        n.reset();
        assert!((n.v_s - (-60.0)).abs() < 1e-10);
        assert!((n.v_d - (-60.0)).abs() < 1e-10);
    }
    #[test]
    fn pr_bounded() {
        let mut n = PinskyRinzelNeuron::new();
        for _ in 0..5000 {
            n.step(50.0, 0.0);
        }
        assert!(n.v_s.is_finite());
    }
    #[test]
    fn pr_dendritic_input() {
        let mut n = PinskyRinzelNeuron::new();
        let _t: i32 = (0..5000).map(|_| n.step(0.0, 5.0)).sum();
        // Dendritic input should also be able to drive spiking
        assert!(n.v_d.is_finite());
    }
    #[test]
    fn pr_nan_no_panic() {
        PinskyRinzelNeuron::new().step(f64::NAN, 0.0);
    }

    // -- HayL5 --
    #[test]
    fn hay_reset() {
        let mut n = HayL5PyramidalNeuron::new();
        for _ in 0..100 {
            n.step(20.0, 0.0);
        }
        n.reset();
        assert!((n.v_s - (-75.0)).abs() < 1e-10);
    }
    #[test]
    fn hay_bounded() {
        let mut n = HayL5PyramidalNeuron::new();
        for _ in 0..500 {
            n.step(100.0, 0.0);
        }
        assert!(n.v_s.is_finite());
    }
    #[test]
    fn hay_nan_no_panic() {
        HayL5PyramidalNeuron::new().step(f64::NAN, 0.0);
    }

    // -- MarderSTG --
    #[test]
    fn marder_reset() {
        let mut n = MarderSTGNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.v - (-60.0)).abs() < 1e-10);
    }
    #[test]
    fn marder_bounded() {
        let mut n = MarderSTGNeuron::new();
        for _ in 0..2000 {
            n.step(50.0);
        }
        assert!(n.v.is_finite());
    }
    #[test]
    fn marder_nan_no_panic() {
        MarderSTGNeuron::new().step(f64::NAN);
    }

    // -- RallCable --
    #[test]
    fn rall_reset() {
        let mut n = RallCableNeuron::new(5);
        for _ in 0..100 {
            n.step(50.0);
        }
        n.reset();
        assert!(n.v.iter().all(|&x| (x - n.v_rest).abs() < 1e-10));
    }
    #[test]
    fn rall_bounded() {
        let mut n = RallCableNeuron::new(5);
        for _ in 0..1000 {
            n.step(500.0);
        }
        assert!(n.v.iter().all(|x| x.is_finite()));
    }
    #[test]
    fn rall_nan_no_panic() {
        RallCableNeuron::new(5).step(f64::NAN);
    }

    // -- BoothRinzel --
    #[test]
    fn booth_reset() {
        let mut n = BoothRinzelNeuron::new();
        for _ in 0..100 {
            n.step(5.0);
        }
        n.reset();
        assert!((n.vs - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn booth_bounded() {
        let mut n = BoothRinzelNeuron::new();
        for _ in 0..2000 {
            n.step(50.0);
        }
        assert!(n.vs.is_finite());
    }
    #[test]
    fn booth_nan_no_panic() {
        BoothRinzelNeuron::new().step(f64::NAN);
    }

    // -- Dendrify --
    #[test]
    fn dendrify_reset() {
        let mut n = DendrifyNeuron::new();
        for _ in 0..100 {
            n.step(50.0);
        }
        n.reset();
        assert!((n.v_s - (-65.0)).abs() < 1e-10);
    }
    #[test]
    fn dendrify_bounded() {
        let mut n = DendrifyNeuron::new();
        for _ in 0..2000 {
            n.step(200.0);
        }
        assert!(n.v_s.is_finite());
    }
    #[test]
    fn dendrify_nan_no_panic() {
        DendrifyNeuron::new().step(f64::NAN);
    }

    // -- TwoCompartmentLIF --
    #[test]
    fn tc_lif_reset() {
        let mut n = TwoCompartmentLIFNeuron::new();
        for _ in 0..50 {
            n.step(0.5, 0.3);
        }
        n.reset();
    }
    #[test]
    fn tc_lif_bounded() {
        let mut n = TwoCompartmentLIFNeuron::new();
        for _ in 0..1000 {
            n.step(100.0, 100.0);
        }
        assert!(n.v_s.is_finite());
    }
    #[test]
    fn tc_lif_nan_no_panic() {
        TwoCompartmentLIFNeuron::new().step(f64::NAN, 0.0);
    }
}

/// Dendritic NMDA spike model.
///
/// Captures the non-linear voltage-dependent Mg²⁺ block of NMDA receptors
/// in dendritic branches. NMDA current has a sigmoidal voltage dependence:
///
///   I_NMDA = g_NMDA · B(V) · (V - E_NMDA)
///   B(V) = 1 / (1 + [Mg²⁺]/3.57 · exp(-0.062 · V))
///
/// This enables coincidence detection: the dendrite only passes current
/// when both presynaptic glutamate AND postsynaptic depolarisation are present.
///
/// Reference: Jahr & Stevens (1990), Schiller et al. (2000).
#[derive(Clone, Debug)]
pub struct DendriticNMDANeuron {
    pub v_soma: f64,
    pub v_dend: f64,
    pub g_nmda: f64,
    pub e_nmda: f64,
    pub mg_conc: f64,
    pub g_coupling: f64,
    pub tau_soma: f64,
    pub tau_dend: f64,
    pub theta: f64,
    pub dt: f64,
}

impl DendriticNMDANeuron {
    pub fn new() -> Self {
        Self {
            v_soma: -65.0,
            v_dend: -65.0,
            g_nmda: 1.5,
            e_nmda: 0.0,
            mg_conc: 1.0,
            g_coupling: 0.5,
            tau_soma: 20.0,
            tau_dend: 50.0,
            theta: -50.0,
            dt: 0.1,
        }
    }

    /// Mg²⁺ block factor (Jahr & Stevens 1990).
    fn mg_block(&self, v: f64) -> f64 {
        1.0 / (1.0 + (self.mg_conc / 3.57) * (-0.062 * v).exp())
    }

    /// Step with somatic input and dendritic glutamate.
    pub fn step(&mut self, i_soma: f64, glutamate: f64) -> i32 {
        // Dendritic NMDA current with Mg²⁺ block.
        let b = self.mg_block(self.v_dend);
        let i_nmda = self.g_nmda * glutamate * b * (self.v_dend - self.e_nmda);

        // Dendrite dynamics.
        let dv_dend =
            (-self.v_dend - 65.0 + i_nmda + self.g_coupling * (self.v_soma - self.v_dend))
                / self.tau_dend;
        self.v_dend += dv_dend * self.dt;

        // Soma dynamics: receives dendritic current + direct input.
        let i_dend_to_soma = self.g_coupling * (self.v_dend - self.v_soma);
        let dv_soma = (-self.v_soma - 65.0 + i_soma + i_dend_to_soma) / self.tau_soma;
        self.v_soma += dv_soma * self.dt;

        if self.v_soma >= self.theta {
            self.v_soma = -65.0;
            1
        } else {
            0
        }
    }

    pub fn reset(&mut self) {
        self.v_soma = -65.0;
        self.v_dend = -65.0;
    }
}

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

/// Multi-compartment neuron (MCN) matching the Spiking-WM architecture.
///
/// Dual-dendrite model with basal and apical compartments. The apical dendrite
/// gates how strongly basal information influences the soma, enabling
/// nonlinear integration for long-term temporal memory in RL tasks.
///
/// Exact equations from arXiv:2503.00713 (Spiking-WM, PNAS 2025):
///
///   τ_b dV_b/dt = -V_b + x_b                                  (basal)
///   τ_a dV_a/dt = -V_a + x_a                                  (apical)
///   τ   dU/dt   = -U + σ(V_a)·[g_B/g_L·(V_b - U) + W_s·I]   (soma)
///   S[t] = Θ(U[t] - V_th)                                     (spike)
///   U[t] ← U[t]·(1 - S[t])                                    (soft reset)
///
/// Default parameters from Table II: τ = τ_a = τ_b = 2.0, g_B/g_L = 1.0,
/// β = 1.0 (sigmoid steepness), V_th = 1.0.
///
/// Reference: Brain-Cog-Lab, arXiv:2503.00713, PNAS 2025.
#[derive(Clone, Debug)]
pub struct MulticompartmentMCNNeuron {
    /// Somatic membrane potential.
    pub u: f64,
    /// Basal dendrite potential.
    pub v_basal: f64,
    /// Apical dendrite potential.
    pub v_apical: f64,
    /// Soma time constant.
    pub tau: f64,
    /// Basal dendrite time constant.
    pub tau_b: f64,
    /// Apical dendrite time constant.
    pub tau_a: f64,
    /// Basal-to-soma conductance ratio (g_B/g_L).
    pub g_ratio: f64,
    /// Sigmoid steepness for apical gating.
    pub beta: f64,
    /// Spike threshold.
    pub v_th: f64,
    /// Time step.
    pub dt: f64,
}

impl MulticompartmentMCNNeuron {
    pub fn new() -> Self {
        Self {
            u: 0.0,
            v_basal: 0.0,
            v_apical: 0.0,
            tau: 2.0,
            tau_b: 2.0,
            tau_a: 2.0,
            g_ratio: 1.0,
            beta: 1.0,
            v_th: 1.0,
            dt: 1.0,
        }
    }

    /// Sigmoid gating function σ(x) = 1/(1 + exp(-βx)).
    fn sigma(&self, x: f64) -> f64 {
        1.0 / (1.0 + (-self.beta * x).exp())
    }

    /// Step with basal input (x_b), apical input (x_a), and direct somatic input.
    pub fn step_compartments(&mut self, x_basal: f64, x_apical: f64, i_soma: f64) -> i32 {
        // Basal dendrite: τ_b dV_b/dt = -V_b + x_b.
        let dv_b = (-self.v_basal + x_basal) / self.tau_b;
        self.v_basal += dv_b * self.dt;

        // Apical dendrite: τ_a dV_a/dt = -V_a + x_a.
        let dv_a = (-self.v_apical + x_apical) / self.tau_a;
        self.v_apical += dv_a * self.dt;

        // Soma: τ dU/dt = -U + σ(V_a)·[g_B/g_L·(V_b - U) + i_soma].
        let gate = self.sigma(self.v_apical);
        let du = (-self.u + gate * (self.g_ratio * (self.v_basal - self.u) + i_soma)) / self.tau;
        self.u += du * self.dt;

        // Spike: S = Θ(U - V_th), reset: U ← U·(1 - S).
        if self.u >= self.v_th {
            self.u = 0.0;
            1
        } else {
            0
        }
    }

    /// Simple step: input goes to basal dendrite only.
    pub fn step(&mut self, current: f64) -> i32 {
        self.step_compartments(current, 0.0, 0.0)
    }

    pub fn reset(&mut self) {
        self.u = 0.0;
        self.v_basal = 0.0;
        self.v_apical = 0.0;
    }
}

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

/// Astrocyte-LIF hybrid unit with calcium wave feedback.
///
/// Models the tripartite synapse: a glial astrocyte monitors extracellular
/// glutamate from a paired LIF neuron and provides slow homeostatic feedback
/// via calcium-dependent gliotransmitter release.
///
///   dCa/dt = -Ca/τ_ca + δ · S_pre(t)        (calcium rise on presynaptic spike)
///   I_glio = g_glio · H(Ca - Ca_thresh)      (gliotransmitter release)
///   dV/dt = -(V - E_L)/τ_m + I_ext + I_glio  (LIF with glial feedback)
///
/// Reference: Perea, Navarrete & Araque, "Tripartite synapses" (2009).
#[derive(Clone, Debug)]
pub struct AstrocyteLIFNeuron {
    pub v: f64,
    pub ca: f64,
    pub tau_m: f64,
    pub tau_ca: f64,
    pub e_l: f64,
    pub theta: f64,
    pub v_reset: f64,
    pub ca_delta: f64,
    pub ca_thresh: f64,
    pub g_glio: f64,
    pub dt: f64,
}

impl AstrocyteLIFNeuron {
    pub fn new() -> Self {
        Self {
            v: -65.0,
            ca: 0.0,
            tau_m: 20.0,
            tau_ca: 500.0,
            e_l: -65.0,
            theta: -50.0,
            v_reset: -65.0,
            ca_delta: 0.1,
            ca_thresh: 0.5,
            g_glio: 2.0,
            dt: 0.1,
        }
    }

    /// Step with external current and presynaptic spike indicator.
    pub fn step_with_pre(&mut self, i_ext: f64, pre_spike: bool) -> i32 {
        // Astrocyte calcium dynamics.
        let dca = -self.ca / self.tau_ca
            + if pre_spike {
                self.ca_delta / self.dt
            } else {
                0.0
            };
        self.ca += dca * self.dt;
        self.ca = self.ca.max(0.0);

        // Gliotransmitter release (Heaviside on calcium).
        let i_glio = if self.ca > self.ca_thresh {
            self.g_glio
        } else {
            0.0
        };

        // LIF membrane dynamics with glial feedback.
        let dv = (-(self.v - self.e_l) + i_ext + i_glio) / self.tau_m;
        self.v += dv * self.dt;

        if self.v >= self.theta {
            self.v = self.v_reset;
            1
        } else {
            0
        }
    }

    /// Simple step (no presynaptic spike).
    pub fn step(&mut self, current: f64) -> i32 {
        self.step_with_pre(current, false)
    }

    pub fn reset(&mut self) {
        self.v = self.e_l;
        self.ca = 0.0;
    }
}

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

// ---- Tests for new multi-compartment / glial models ----

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

    #[test]
    fn nmda_coincidence_detection() {
        let mut n = DendriticNMDANeuron::new();
        // Only soma input — dendrite contributes little.
        let mut spikes_soma_only = 0;
        for _ in 0..2000 {
            spikes_soma_only += n.step(8.0, 0.0);
        }
        n.reset();
        // Soma + glutamate — NMDA amplifies.
        let mut spikes_both = 0;
        for _ in 0..2000 {
            spikes_both += n.step(8.0, 1.0);
        }
        // Coincidence: both inputs together should fire more.
        assert!(
            spikes_both >= spikes_soma_only,
            "NMDA coincidence: both={spikes_both} must >= soma_only={spikes_soma_only}"
        );
    }

    #[test]
    fn nmda_mg_block_voltage_dependent() {
        let n = DendriticNMDANeuron::new();
        let b_hyper = n.mg_block(-80.0);
        let b_depol = n.mg_block(-20.0);
        assert!(
            b_depol > b_hyper,
            "Mg block must relieve at depolarised potentials: B(-20)={b_depol:.3} > B(-80)={b_hyper:.3}"
        );
    }

    #[test]
    fn nmda_zero_glutamate_no_nmda_current() {
        let mut n = DendriticNMDANeuron::new();
        let spikes: i32 = (0..500).map(|_| n.step(0.0, 0.0)).sum();
        assert_eq!(spikes, 0, "No input → no spikes");
    }

    #[test]
    fn mcn_apical_gating() {
        // Without apical input, gate = σ(0) = 0.5, moderate drive.
        let mut n_no_apical = MulticompartmentMCNNeuron::new();
        let mut spikes_no = 0;
        for _ in 0..100 {
            spikes_no += n_no_apical.step_compartments(3.0, 0.0, 0.0);
        }
        // With strong apical input, gate ≈ 1.0, full basal→soma coupling.
        let mut n_apical = MulticompartmentMCNNeuron::new();
        let mut spikes_yes = 0;
        for _ in 0..100 {
            spikes_yes += n_apical.step_compartments(3.0, 5.0, 0.0);
        }
        assert!(
            spikes_yes >= spikes_no,
            "Apical gating should boost firing: apical={spikes_yes} >= none={spikes_no}"
        );
    }

    #[test]
    fn mcn_basal_dendrite_memory() {
        // τ_b = 2.0, dt = 1.0: V_b decays by factor (1 - dt/τ) = 0.5 per step.
        let mut n = MulticompartmentMCNNeuron::new();
        n.step_compartments(5.0, 0.0, 0.0);
        let v_after = n.v_basal;
        n.step_compartments(0.0, 0.0, 0.0);
        let v_decay = n.v_basal;
        assert!(
            v_decay.abs() > 0.1 * v_after.abs(),
            "Basal dendrite retains memory: {v_decay:.3} vs {v_after:.3}"
        );
    }

    #[test]
    fn mcn_reset_clears_all() {
        let mut n = MulticompartmentMCNNeuron::new();
        for _ in 0..50 {
            n.step(2.0);
        }
        n.reset();
        assert_eq!(n.u, 0.0);
        assert_eq!(n.v_basal, 0.0);
        assert_eq!(n.v_apical, 0.0);
    }

    #[test]
    fn astrocyte_calcium_rises_on_pre_spikes() {
        let mut n = AstrocyteLIFNeuron::new();
        let ca_before = n.ca;
        for _ in 0..100 {
            n.step_with_pre(0.0, true);
        }
        assert!(
            n.ca > ca_before,
            "Calcium must rise with presynaptic spikes"
        );
    }

    #[test]
    fn astrocyte_gliotransmitter_boosts_firing() {
        let mut n_no_glio = AstrocyteLIFNeuron::new();
        let mut n_glio = AstrocyteLIFNeuron::new();

        let mut spikes_no = 0;
        let mut spikes_yes = 0;
        for _ in 0..5000 {
            spikes_no += n_no_glio.step_with_pre(10.0, false);
            spikes_yes += n_glio.step_with_pre(10.0, true); // pre spikes → Ca → glio
        }
        assert!(
            spikes_yes >= spikes_no,
            "Gliotransmitter should boost firing: with={spikes_yes} >= without={spikes_no}"
        );
    }

    #[test]
    fn astrocyte_calcium_decays() {
        let mut n = AstrocyteLIFNeuron::new();
        // Build up calcium.
        for _ in 0..200 {
            n.step_with_pre(0.0, true);
        }
        let ca_peak = n.ca;
        // Let it decay.
        for _ in 0..5000 {
            n.step_with_pre(0.0, false);
        }
        assert!(
            n.ca < ca_peak * 0.5,
            "Calcium must decay: current={:.4} < peak={:.4}*0.5",
            n.ca,
            ca_peak
        );
    }
}