sc_neurocore_engine/neurons/interneurons/
vip_neuron.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later
2// Commercial license available
3// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
4// © Code 2020–2026 Miroslav Šotek. All rights reserved.
5// ORCID: 0009-0009-3560-0851
6// Contact: www.anulum.li | protoscience@anulum.li
7// SC-NeuroCore — VIP interneuron model
8
9/// VIP (vasoactive intestinal peptide) irregular-spiking interneuron.
10///
11/// Biophysics: Na+, K+, A-type K+ (Kv4, transient outward, causes
12/// accommodation), leak. High input resistance, small soma.
13/// Key properties: irregular/accommodating firing, disinhibitory
14/// role (inhibits SST+ and PV+), bipolar morphology.
15///
16/// Based on Porter et al. 1998 / Bhatt et al. 2019 parameterisation.
17#[derive(Clone, Debug)]
18pub struct VIPNeuron {
19    pub v: f64,
20    pub h: f64,
21    pub n: f64,
22    pub a: f64, // A-type K+ activation
23    pub b: f64, // A-type K+ inactivation
24    // Conductances
25    pub g_na: f64,
26    pub g_k: f64,
27    pub g_a: f64,
28    pub g_l: f64,
29    // Reversal potentials
30    pub e_na: f64,
31    pub e_k: f64,
32    pub e_l: f64,
33    pub c_m: f64,
34    pub dt: f64,
35    pub v_threshold: f64,
36}
37
38impl VIPNeuron {
39    pub fn new() -> Self {
40        Self {
41            v: -65.0,
42            h: 0.8,
43            n: 0.1,
44            a: 0.0,
45            b: 0.9,
46            g_na: 35.0, // Lower than PV+ (smaller soma)
47            g_k: 6.0,
48            g_a: 8.0,  // Strong A-current → accommodation
49            g_l: 0.01, // High input resistance
50            e_na: 55.0,
51            e_k: -90.0,
52            e_l: -65.0,
53            c_m: 0.5, // Small soma → low capacitance
54            dt: 0.025,
55            v_threshold: -20.0,
56        }
57    }
58
59    /// Return `[dV, dh, dn, da, db]` of the five-state VIP system at one consistent
60    /// state. All gates relax through sigmoidal steady states (no singularities).
61    fn derivatives(&self, v: f64, h: f64, n: f64, a: f64, b: f64, current: f64) -> [f64; 5] {
62        let m_inf = 1.0 / (1.0 + (-(v + 30.0) / 9.5).exp());
63        let h_inf = 1.0 / (1.0 + ((v + 53.0) / 7.0).exp());
64        let tau_h = 0.37 + 2.78 / (1.0 + ((v + 40.5) / 6.0).exp());
65        let n_inf = 1.0 / (1.0 + (-(v + 30.0) / 10.0).exp());
66        let tau_n = 0.37 + 1.85 / (1.0 + ((v + 27.0) / 15.0).exp());
67        let a_inf = 1.0 / (1.0 + (-(v + 50.0) / 20.0).exp());
68        let b_inf = 1.0 / (1.0 + ((v + 78.0) / 6.0).exp());
69        let dh = (h_inf - h) / tau_h;
70        let dn = (n_inf - n) / tau_n;
71        let da = (a_inf - a) / 5.0;
72        let db = (b_inf - b) / 50.0;
73        let i_na = self.g_na * m_inf * m_inf * m_inf * h * (v - self.e_na);
74        let i_k = self.g_k * n * n * n * n * (v - self.e_k);
75        let i_a = self.g_a * a * a * a * b * (v - self.e_k);
76        let i_l = self.g_l * (v - self.e_l);
77        let dv = (-i_na - i_k - i_a - i_l + current) / self.c_m;
78        [dv, dh, dn, da, db]
79    }
80
81    /// Return one classical RK4 increment of `[V, h, n, a, b]`, holding `current`
82    /// constant across the four stages.
83    fn rk4_substep(&self, s: [f64; 5], current: f64) -> [f64; 5] {
84        let dt = self.dt;
85        let k1 = self.derivatives(s[0], s[1], s[2], s[3], s[4], current);
86        let k2 = self.derivatives(
87            s[0] + 0.5 * dt * k1[0],
88            s[1] + 0.5 * dt * k1[1],
89            s[2] + 0.5 * dt * k1[2],
90            s[3] + 0.5 * dt * k1[3],
91            s[4] + 0.5 * dt * k1[4],
92            current,
93        );
94        let k3 = self.derivatives(
95            s[0] + 0.5 * dt * k2[0],
96            s[1] + 0.5 * dt * k2[1],
97            s[2] + 0.5 * dt * k2[2],
98            s[3] + 0.5 * dt * k2[3],
99            s[4] + 0.5 * dt * k2[4],
100            current,
101        );
102        let k4 = self.derivatives(
103            s[0] + dt * k3[0],
104            s[1] + dt * k3[1],
105            s[2] + dt * k3[2],
106            s[3] + dt * k3[3],
107            s[4] + dt * k3[4],
108            current,
109        );
110        let mut out = [0.0_f64; 5];
111        for i in 0..5 {
112            out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
113        }
114        out
115    }
116
117    pub fn step(&mut self, current: f64) -> i32 {
118        let v_prev = self.v;
119        let mut s = [self.v, self.h, self.n, self.a, self.b];
120        for _ in 0..4 {
121            s = self.rk4_substep(s, current);
122        }
123        self.v = s[0];
124        self.h = s[1];
125        self.n = s[2];
126        self.a = s[3];
127        self.b = s[4];
128        if self.v >= self.v_threshold && v_prev < self.v_threshold {
129            1
130        } else {
131            0
132        }
133    }
134
135    pub fn reset(&mut self) {
136        self.v = -65.0;
137        self.h = 0.8;
138        self.n = 0.1;
139        self.a = 0.0;
140        self.b = 0.9;
141    }
142}
143
144impl Default for VIPNeuron {
145    fn default() -> Self {
146        Self::new()
147    }
148}
149
150// ═══════════════════════════════════════════════════════════════════
151// Chandelier Cell (Axo-Axonic)
152// ═══════════════════════════════════════════════════════════════════
153
154#[cfg(test)]
155mod tests {
156    use super::*;
157
158    #[test]
159    fn vip_fires_with_input() {
160        let mut n = VIPNeuron::new();
161        let spikes: i32 = (0..10000).map(|_| n.step(2.0)).sum();
162        assert!(spikes > 0, "VIP must fire with sustained input");
163    }
164
165    #[test]
166    fn vip_no_fire_without_input() {
167        let mut n = VIPNeuron::new();
168        let spikes: i32 = (0..5000).map(|_| n.step(0.0)).sum();
169        assert_eq!(spikes, 0);
170    }
171
172    #[test]
173    fn vip_accommodation() {
174        // A-current causes transient accommodation at spike onset.
175        // Compare fresh neuron's first 100 steps vs steady-state.
176        let mut n = VIPNeuron::new();
177        // First 500 steps: A-current b gate is high → strong IA → suppresses early spikes
178        let onset: i32 = (0..500).map(|_| n.step(3.0)).sum();
179        // Skip 5000 steps to reach steady state
180        for _ in 0..5000 {
181            n.step(3.0);
182        }
183        // Next 500 steps at steady state
184        let steady: i32 = (0..500).map(|_| n.step(3.0)).sum();
185        // At steady state, b has dropped, IA is weaker → fires at least as much
186        assert!(
187            steady >= onset,
188            "VIP steady-state ({steady}) should fire >= onset ({onset})"
189        );
190    }
191
192    #[test]
193    fn vip_reset_roundtrip() {
194        let mut n = VIPNeuron::new();
195        for _ in 0..5000 {
196            n.step(3.0);
197        }
198        n.reset();
199        let mut fresh = VIPNeuron::new();
200        let r1: i32 = (0..2000).map(|_| n.step(3.0)).sum();
201        let r2: i32 = (0..2000).map(|_| fresh.step(3.0)).sum();
202        assert_eq!(r1, r2);
203    }
204
205    #[test]
206    fn vip_voltage_bounded() {
207        let mut n = VIPNeuron::new();
208        for _ in 0..20000 {
209            n.step(5.0);
210        }
211        assert!(n.v.is_finite());
212    }
213
214    #[test]
215    #[ignore = "wall-clock performance smoke; use Criterion benches for timing evidence"]
216    fn vip_performance_10k_steps() {
217        let mut n = VIPNeuron::new();
218        let start = std::time::Instant::now();
219        for _ in 0..10_000 {
220            n.step(3.0);
221        }
222        assert!(start.elapsed().as_millis() < 100);
223    }
224}
sc_neurocore_engine/neurons/interneurons/vip_neuron.rs

sc_neurocore_engine/neurons/interneurons/
vip_neuron.rs
