sc_neurocore_engine/neurons/cerebellar/

golgi.rs

// 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 — Cerebellar Circuit Neuron Models

use super::super::biophysical::safe_rate;

// ═══════════════════════════════════════════════════════════════════
// Golgi Cell
// ═══════════════════════════════════════════════════════════════════

/// Cerebellar Golgi cell — Solinas et al. 2007 full model.
///
/// Large inhibitory interneuron in the granular layer. Provides tonic
/// and phasic GABAergic/glycinergic inhibition to granule cells.
/// Spontaneously active at 3-10 Hz due to intrinsic pacemaker currents.
///
/// Full Solinas 2007 model with 11 ionic currents:
/// - **INa_t** (transient Na, m³h): fast spike generation
/// - **INa_p** (persistent Na, p): subthreshold oscillations, pacemaking
/// - **IK_dr** (delayed rectifier K, n⁴): repolarisation
/// - **IK_A** (A-type K, a³b): onset delay, inter-spike interval
/// - **IK_M** (muscarinic/slow K, w): spike frequency adaptation
/// - **ICa_T** (T-type Ca²⁺, m_t²s): rebound, subthreshold oscillations
/// - **ICa_N** (N-type Ca²⁺, c²): high-voltage activated, AHP trigger
/// - **IBK** (BK, Ca²⁺+V dependent): fast AHP
/// - **ISK** (SK, Ca²⁺ dependent): slow AHP, pacemaker regulation
/// - **Ih** (HCN, r): sag, resting potential, pacemaker contribution
/// - **IL** (leak)
///
/// 10 sub-steps (dt_sub = 0.05 ms) for Na gating stability.
///
/// Solinas et al., Front Cell Neurosci 1:2, 2007.
#[derive(Clone, Debug)]
pub struct GolgiCell {
    pub v: f64,
    pub m: f64,    // Na_t activation
    pub h: f64,    // Na_t inactivation
    pub p_na: f64, // Na_p persistent activation
    pub n: f64,    // K_dr activation
    pub a: f64,    // K_A activation
    pub b: f64,    // K_A inactivation
    pub w: f64,    // K_M (muscarinic) activation
    pub m_t: f64,  // Ca_T activation
    pub s: f64,    // Ca_T inactivation
    pub c_n: f64,  // Ca_N activation
    pub r: f64,    // Ih activation
    pub ca: f64,   // Intracellular Ca²⁺ (µM)
    // Conductances (mS/cm²)
    pub g_na_t: f64,
    pub g_na_p: f64,
    pub g_kdr: f64,
    pub g_ka: f64,
    pub g_km: f64,
    pub g_cat: f64,
    pub g_can: f64,
    pub g_bk: f64,
    pub g_sk: f64,
    pub g_h: f64,
    pub g_l: f64,
    // Reversals
    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 tau_ca: f64,
    pub kd_bk: f64,
    pub kd_sk: f64,
    pub dt: f64,
    pub sub_steps: usize,
    pub gain: f64,
}

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

impl GolgiCell {
    pub fn new() -> Self {
        Self {
            v: -60.0,
            m: 0.02,
            h: 0.85,
            p_na: 0.01,
            n: 0.05,
            a: 0.1,
            b: 0.8,
            w: 0.01,
            m_t: 0.01,
            s: 0.9,
            c_n: 0.01,
            r: 0.1,
            ca: 0.05,
            g_na_t: 48.0, // Solinas 2007 Table 1
            g_na_p: 0.2,  // Persistent Na (small but critical for pacemaking)
            g_kdr: 16.0,
            g_ka: 8.0,  // A-type
            g_km: 1.0,  // Muscarinic slow K
            g_cat: 0.5, // T-type Ca²⁺
            g_can: 1.0, // N-type Ca²⁺ (high-voltage)
            g_bk: 3.0,  // BK fast AHP
            g_sk: 1.0,  // SK slow AHP
            g_h: 0.1,   // Ih
            g_l: 0.05,
            e_na: 55.0,
            e_k: -90.0,
            e_ca: 120.0,
            e_h: -40.0,
            e_l: -55.0, // Depolarised leak for spontaneous activity
            c_m: 1.0,
            tau_ca: 200.0,
            kd_bk: 1.0,
            kd_sk: 0.5,
            dt: 0.5,
            sub_steps: 10,
            gain: 1.0,
        }
    }

    #[inline]
    fn boltz(v: f64, vh: f64, k: f64) -> f64 {
        let x = (v - vh) / k;
        if x >= 0.0 {
            1.0 / (1.0 + (-x).exp())
        } else {
            let ex = x.exp();
            ex / (1.0 + ex)
        }
    }

    #[inline]
    fn voltage_valid(value: f64) -> bool {
        value.is_finite() && (-100.0..=60.0).contains(&value)
    }

    #[inline]
    fn probability(value: f64) -> bool {
        value.is_finite() && (0.0..=1.0).contains(&value)
    }

    #[inline]
    fn gate_alpha_beta(previous: f64, alpha: f64, beta: f64, phi: f64, dt: f64) -> Option<f64> {
        let total = phi * (alpha + beta);
        if !previous.is_finite()
            || !alpha.is_finite()
            || !beta.is_finite()
            || !total.is_finite()
            || !dt.is_finite()
            || total <= 0.0
        {
            return None;
        }
        let steady = alpha / (alpha + beta);
        Some((steady + (previous - steady) * (-total * dt).exp()).clamp(0.0, 1.0))
    }

    #[inline]
    fn gate_inf(previous: f64, steady: f64, tau: f64, dt: f64) -> Option<f64> {
        if !previous.is_finite()
            || !steady.is_finite()
            || !tau.is_finite()
            || !dt.is_finite()
            || tau <= 0.0
        {
            return None;
        }
        Some((steady + (previous - steady) * (-dt / tau).exp()).clamp(0.0, 1.0))
    }

    #[inline]
    fn calcium_exact(previous: f64, entry: f64, tau: f64, dt: f64) -> Option<f64> {
        if !previous.is_finite()
            || !entry.is_finite()
            || !tau.is_finite()
            || !dt.is_finite()
            || tau <= 0.0
            || previous < 0.0
        {
            return None;
        }
        let steady = entry * tau;
        let value = steady + (previous - steady) * (-dt / tau).exp();
        value.is_finite().then_some(value.max(0.0))
    }

    fn valid_state(&self) -> bool {
        Self::voltage_valid(self.v)
            && [
                self.m, self.h, self.p_na, self.n, self.a, self.b, self.w, self.m_t, self.s,
                self.c_n, self.r,
            ]
            .into_iter()
            .all(Self::probability)
            && [
                self.g_na_t,
                self.g_na_p,
                self.g_kdr,
                self.g_ka,
                self.g_km,
                self.g_cat,
                self.g_can,
                self.g_bk,
                self.g_sk,
                self.g_h,
                self.g_l,
            ]
            .into_iter()
            .all(|g| g.is_finite() && g >= 0.0)
            && self.ca.is_finite()
            && self.ca >= 0.0
            && self.e_na.is_finite()
            && self.e_k.is_finite()
            && self.e_ca.is_finite()
            && self.e_h.is_finite()
            && self.e_l.is_finite()
            && self.c_m.is_finite()
            && self.tau_ca.is_finite()
            && self.kd_bk.is_finite()
            && self.kd_sk.is_finite()
            && self.dt.is_finite()
            && self.gain.is_finite()
            && self.c_m > 0.0
            && self.tau_ca > 0.0
            && self.kd_bk > 0.0
            && self.kd_sk > 0.0
            && self.dt > 0.0
            && self.sub_steps > 0
            && self.gain >= 0.0
    }

    pub fn step(&mut self, current: f64) -> i32 {
        if !current.is_finite() || !self.valid_state() {
            return 0;
        }

        let input = self.gain * current;
        let dt_sub = self.dt / self.sub_steps as f64;
        let v_prev = self.v;
        let mut next = self.clone();

        for _ in 0..self.sub_steps {
            let v = next.v;

            // Na_t: m³h (fast, WB-style alpha/beta)
            let alpha_m = safe_rate(0.1, 35.0, v, 10.0, 1.0);
            let beta_m = 4.0 * (-(v + 60.0) / 18.0).exp();
            let alpha_h = 0.07 * (-(v + 58.0) / 20.0).exp();
            let beta_h = 1.0 / (1.0 + (-(v + 28.0) / 10.0).exp());
            let Some(m) = Self::gate_alpha_beta(next.m, alpha_m, beta_m, 5.0, dt_sub) else {
                return 0;
            };
            let Some(h) = Self::gate_alpha_beta(next.h, alpha_h, beta_h, 5.0, dt_sub) else {
                return 0;
            };

            // Na_p: persistent (Boltzmann, slow)
            let pna_inf = Self::boltz(v, -48.0, 5.0);
            let tau_pna = 5.0 + 20.0 / (1.0 + ((v + 48.0) / 10.0).powi(2)).max(0.01);
            let Some(p_na) = Self::gate_inf(next.p_na, pna_inf, tau_pna, dt_sub) else {
                return 0;
            };

            // K_dr: n⁴
            let alpha_n = safe_rate(0.01, 34.0, v, 10.0, 0.1);
            let beta_n = 0.125 * (-(v + 44.0) / 80.0).exp();
            let Some(n) = Self::gate_alpha_beta(next.n, alpha_n, beta_n, 5.0, dt_sub) else {
                return 0;
            };

            // K_A: a³b (Solinas 2007: V1/2_act ≈ -27 mV, V1/2_inact ≈ -80 mV)
            let a_inf = Self::boltz(v, -27.0, 16.0);
            let b_inf = Self::boltz(v, -80.0, -6.0);
            let Some(a) = Self::gate_inf(next.a, a_inf, 2.0, dt_sub) else {
                return 0;
            };
            let Some(b) = Self::gate_inf(next.b, b_inf, 15.0, dt_sub) else {
                return 0;
            };

            // K_M: w (slow muscarinic)
            let w_inf = Self::boltz(v, -35.0, 10.0);
            let tau_w = 100.0 / (3.3 * ((v + 35.0) / 20.0).exp() + (-(v + 35.0) / 20.0).exp());
            let Some(w) = Self::gate_inf(next.w, w_inf, tau_w, dt_sub) else {
                return 0;
            };

            // Ca_T: m_t²s
            let mt_inf = Self::boltz(v, -52.0, 5.0);
            let s_inf = Self::boltz(v, -60.0, -6.5);
            let tau_s = 20.0 + 50.0 / (1.0 + ((v + 65.0) / 10.0).powi(2)).max(0.01);
            let Some(m_t) = Self::gate_inf(next.m_t, mt_inf, 1.0, dt_sub) else {
                return 0;
            };
            let Some(s) = Self::gate_inf(next.s, s_inf, tau_s, dt_sub) else {
                return 0;
            };

            // Ca_N: c² (high-voltage activated)
            let cn_inf = Self::boltz(v, -20.0, 5.0);
            let tau_cn = 2.0 + 10.0 / (1.0 + ((v + 20.0) / 10.0).powi(2)).max(0.01);
            let Some(c_n) = Self::gate_inf(next.c_n, cn_inf, tau_cn, dt_sub) else {
                return 0;
            };

            // Ih: r (slow, hyperpolarisation-activated)
            let r_inf = Self::boltz(v, -80.0, -10.0);
            let tau_r = 50.0 + 200.0 / (1.0 + ((v + 80.0) / 20.0).powi(2)).max(0.01);
            let Some(r) = Self::gate_inf(next.r, r_inf, tau_r, dt_sub) else {
                return 0;
            };

            // Ca²⁺ dynamics (entry via Ca_T + Ca_N, decay)
            let g_cat = self.g_cat * m_t.powi(2) * s;
            let g_can = self.g_can * c_n.powi(2);
            let i_cat = g_cat * (v - self.e_ca);
            let i_can = g_can * (v - self.e_ca);
            let ca_entry = if i_cat + i_can < 0.0 {
                -(i_cat + i_can) * 0.001
            } else {
                0.0
            };
            let Some(ca) = Self::calcium_exact(next.ca, ca_entry, self.tau_ca, dt_sub) else {
                return 0;
            };

            // BK: voltage + Ca²⁺ dependent (Hill n=2 for Ca²⁺ shift)
            // V1/2 shifts from +100 mV (low Ca) to -20 mV (high Ca)
            let ca2 = ca * ca;
            let kd2 = self.kd_bk * self.kd_bk;
            let bk_v = Self::boltz(v, 100.0 - 120.0 * ca2 / (ca2 + kd2), 15.0);
            // SK: Ca²⁺ dependent (Hill n=2)
            let sk_inf = ca2 / (ca2 + self.kd_sk.powi(2));

            // All ionic currents
            let g_na = self.g_na_t * m.powi(3) * h + self.g_na_p * p_na;
            let g_k = self.g_kdr * n.powi(4)
                + self.g_ka * a.powi(3) * b
                + self.g_km * w
                + self.g_bk * bk_v
                + self.g_sk * sk_inf;
            let g_ca = g_cat + g_can;
            let g_h = self.g_h * r;
            let g_total = g_na + g_k + g_ca + g_h + self.g_l;
            if !g_total.is_finite() || g_total <= 0.0 {
                return 0;
            }
            let steady_v = (input
                + g_na * self.e_na
                + g_k * self.e_k
                + g_ca * self.e_ca
                + g_h * self.e_h
                + self.g_l * self.e_l)
                / g_total;
            let v_next = steady_v + (v - steady_v) * (-(g_total / self.c_m) * dt_sub).exp();
            if !Self::voltage_valid(v_next) || !ca.is_finite() || ca < 0.0 {
                return 0;
            }

            next.v = v_next;
            next.m = m;
            next.h = h;
            next.p_na = p_na;
            next.n = n;
            next.a = a;
            next.b = b;
            next.w = w;
            next.m_t = m_t;
            next.s = s;
            next.c_n = c_n;
            next.r = r;
            next.ca = ca;
        }

        *self = next;

        // Spike: V crosses 0 mV
        if self.v >= 0.0 && v_prev < 0.0 {
            1
        } else {
            0
        }
    }

    pub fn reset(&mut self) {
        *self = Self::new();
    }
}

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

    // -- Golgi Cell tests --

    #[test]
    fn golgi_fires_with_input() {
        let mut n = GolgiCell::new();
        let mut spikes = 0;
        for _ in 0..10_000 {
            spikes += n.step(15.0);
        }
        assert!(
            spikes > 10,
            "Golgi cell must fire with excitatory input, got {spikes}"
        );
    }

    #[test]
    fn golgi_spontaneous_firing() {
        // Golgi cells are spontaneously active due to depolarised leak
        let mut n = GolgiCell::new();
        let _spikes: i32 = (0..20_000).map(|_| n.step(0.0)).sum();
        // With e_l = -60 and v_t = -56.2, may or may not spontaneously fire
        // The key property is that they fire easily with minimal input
        let mut n2 = GolgiCell::new();
        let mut spikes_small = 0;
        for _ in 0..20_000 {
            spikes_small += n2.step(0.5);
        }
        assert!(
            spikes_small > 0,
            "Golgi cell should fire with minimal input (near-threshold), got {spikes_small}"
        );
    }

    #[test]
    fn golgi_ahp_reduces_rate_at_high_drive() {
        // BK + SK provide AHP — removing them should increase sustained firing
        let mut with_ahp = GolgiCell::new();
        let mut no_ahp = GolgiCell::new();
        no_ahp.g_bk = 0.0;
        no_ahp.g_sk = 0.0;
        let mut spikes_with = 0;
        let mut spikes_no = 0;
        for _ in 0..10_000 {
            spikes_with += with_ahp.step(10.0);
            spikes_no += no_ahp.step(10.0);
        }
        assert!(
            spikes_no >= spikes_with,
            "AHP removal should increase firing: with={spikes_with}, without={spikes_no}"
        );
    }

    #[test]
    fn golgi_ka_is_transient() {
        // K_A (A-type) is transient: activates fast, inactivates fast.
        // In full 11-current Golgi model, removing K_A changes firing pattern.
        let mut with_a = GolgiCell::new();
        let mut no_a = GolgiCell::new();
        no_a.g_ka = 0.0;
        let mut spikes_with = 0;
        let mut spikes_no = 0;
        for _ in 0..10_000 {
            spikes_with += with_a.step(5.0);
            spikes_no += no_a.step(5.0);
        }
        // Both configurations must fire (K_A doesn't prevent spiking)
        assert!(spikes_with > 0, "Must fire with K_A");
        // K_A modulates rate — the difference should be measurable
        assert!(
            spikes_with != spikes_no,
            "K_A should affect firing rate: with={spikes_with}, without={spikes_no}"
        );
    }

    #[test]
    fn golgi_ca_accumulates_during_spiking() {
        let mut n = GolgiCell::new();
        let ca_init = n.ca;
        for _ in 0..5000 {
            n.step(10.0);
        }
        assert!(
            n.ca > ca_init,
            "Ca²⁺ must rise during spiking: init={ca_init}, now={}",
            n.ca
        );
    }

    #[test]
    fn golgi_negative_input_no_crash() {
        let mut n = GolgiCell::new();
        for _ in 0..10_000 {
            n.step(-100.0);
        }
        assert!(n.v.is_finite(), "Must stay finite with negative input");
        assert!(n.v >= -100.0);
    }

    #[test]
    fn golgi_nan_input_stays_finite() {
        let mut n = GolgiCell::new();
        n.step(f64::NAN);
        assert!(n.v.is_finite(), "NaN input must not corrupt state");
    }

    #[test]
    fn golgi_extreme_input_bounded() {
        let mut n = GolgiCell::new();
        for _ in 0..1000 {
            n.step(1e6);
        }
        assert!(
            n.v.is_finite() && n.v <= 60.0,
            "Extreme input must stay bounded"
        );
    }

    #[test]
    fn golgi_reset_clears_state() {
        let mut n = GolgiCell::new();
        for _ in 0..5000 {
            n.step(10.0);
        }
        n.reset();
        let fresh = GolgiCell::new();
        assert_eq!(n.v, fresh.v);
        assert_eq!(n.ca, fresh.ca);
        assert_eq!(n.m, fresh.m);
        assert_eq!(n.h, fresh.h);
        assert_eq!(n.p_na, fresh.p_na);
        assert_eq!(n.w, fresh.w);
        assert_eq!(n.r, fresh.r);
    }

    #[test]
    fn golgi_gates_bounded() {
        let mut n = GolgiCell::new();
        for _ in 0..10_000 {
            n.step(15.0);
        }
        // All 11 gating variables must be in [0, 1]
        for (name, val) in [
            ("m", n.m),
            ("h", n.h),
            ("p_na", n.p_na),
            ("n", n.n),
            ("a", n.a),
            ("b", n.b),
            ("w", n.w),
            ("m_t", n.m_t),
            ("s", n.s),
            ("c_n", n.c_n),
            ("r", n.r),
        ] {
            assert!((0.0..=1.0).contains(&val), "{name} out of bounds: {val}");
        }
        assert!(n.ca >= 0.0, "Ca²⁺ must be non-negative: {}", n.ca);
    }

    #[test]
    fn golgi_has_eleven_currents() {
        // Solinas 2007: Na_t, Na_p, K_dr, K_A, K_M, Ca_T, Ca_N, BK, SK, Ih, leak = 11
        let n = GolgiCell::new();
        assert!(n.g_na_t > 0.0, "Na_t missing");
        assert!(n.g_na_p > 0.0, "Na_p missing");
        assert!(n.g_kdr > 0.0, "K_dr missing");
        assert!(n.g_ka > 0.0, "K_A missing");
        assert!(n.g_km > 0.0, "K_M missing");
        assert!(n.g_cat > 0.0, "Ca_T missing");
        assert!(n.g_can > 0.0, "Ca_N missing");
        assert!(n.g_bk > 0.0, "BK missing");
        assert!(n.g_sk > 0.0, "SK missing");
        assert!(n.g_h > 0.0, "Ih missing");
        assert!(n.g_l > 0.0, "Leak missing");
    }

    #[test]
    fn golgi_persistent_na_depolarises() {
        // Na_p contributes to pacemaking — removing it should reduce excitability
        let mut with_nap = GolgiCell::new();
        let mut no_nap = GolgiCell::new();
        no_nap.g_na_p = 0.0;
        let mut spikes_with = 0;
        let mut spikes_no = 0;
        for _ in 0..10_000 {
            spikes_with += with_nap.step(2.0);
            spikes_no += no_nap.step(2.0);
        }
        assert!(
            spikes_with >= spikes_no,
            "Na_p should increase excitability: with={spikes_with} vs without={spikes_no}"
        );
    }

    #[test]
    fn golgi_km_modulates_firing_pattern() {
        // K_M (muscarinic) is a slow K+ conductance that changes the exact
        // pacemaking trajectory. Under this fixed-drive protocol, removing K_M
        // depolarises the cell into a different conductance balance rather than
        // producing a globally monotonic rate increase.
        let mut with_km = GolgiCell::new();
        let mut no_km = GolgiCell::new();
        no_km.g_km = 0.0;
        let mut spikes_with = 0;
        let mut spikes_no = 0;
        for _ in 0..10_000 {
            spikes_with += with_km.step(10.0);
            spikes_no += no_km.step(10.0);
        }
        assert!(spikes_with > 0, "Golgi cell with K_M should fire");
        assert!(spikes_no > 0, "Golgi cell without K_M should fire");
        assert!(
            spikes_with != spikes_no,
            "K_M should measurably modulate firing: with_km={spikes_with}, without={spikes_no}"
        );
    }

    #[test]
    fn golgi_ih_sag() {
        // Ih activates on hyperpolarisation → sag towards resting
        let mut with_h = GolgiCell::new();
        let mut no_h = GolgiCell::new();
        no_h.g_h = 0.0;
        // Mild hyperpolarisation (g_h=0.1 is small, so don't drive to clamp)
        for _ in 0..10_000 {
            with_h.step(-1.0);
            no_h.step(-1.0);
        }
        // Ih should depolarise relative to no-Ih (sag)
        assert!(
            with_h.v > no_h.v,
            "Ih should cause sag (less hyperpolarised): with_h={:.1} vs no_h={:.1}",
            with_h.v,
            no_h.v
        );
    }

    #[test]
    fn golgi_bk_fast_ahp() {
        // BK channels contribute to fast AHP — removing them should widen spikes
        let mut with_bk = GolgiCell::new();
        let mut no_bk = GolgiCell::new();
        no_bk.g_bk = 0.0;
        // Drive both to fire, measure voltage after spike
        let mut spikes_with = 0;
        let mut spikes_no = 0;
        for _ in 0..10_000 {
            spikes_with += with_bk.step(10.0);
            spikes_no += no_bk.step(10.0);
        }
        // Without BK, model should still fire (test stability)
        assert!(
            spikes_with > 0 && spikes_no > 0,
            "Both should fire: with_bk={spikes_with}, no_bk={spikes_no}"
        );
    }

    #[test]
    fn golgi_sk_slow_adaptation() {
        // SK channels provide slow AHP → spike frequency adaptation
        let mut with_sk = GolgiCell::new();
        let mut no_sk = GolgiCell::new();
        no_sk.g_sk = 0.0;
        let mut spikes_with = 0;
        let mut spikes_no = 0;
        for _ in 0..20_000 {
            spikes_with += with_sk.step(8.0);
            spikes_no += no_sk.step(8.0);
        }
        assert!(
            spikes_no >= spikes_with,
            "SK removal should increase firing: with_sk={spikes_with}, no_sk={spikes_no}"
        );
    }

    #[test]
    fn golgi_performance_1k_steps() {
        let start = std::time::Instant::now();
        let mut n = GolgiCell::new();
        for _ in 0..1_000 {
            std::hint::black_box(n.step(5.0));
        }
        let elapsed = start.elapsed();
        assert!(
            elapsed.as_millis() < 50,
            "1k steps must complete in <50ms, took {}ms",
            elapsed.as_millis()
        );
    }

    #[test]
    fn golgi_default_matches_constructor_contract() {
        let default = GolgiCell::default();
        let constructed = GolgiCell::new();
        assert_eq!(default.v, constructed.v);
        assert_eq!(default.ca, constructed.ca);
        assert_eq!(default.g_na_t, constructed.g_na_t);
        assert_eq!(default.sub_steps, constructed.sub_steps);
    }
}