sc_neurocore_engine/neurons/interneurons/
martinotti_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 — Martinotti interneuron model
8
9/// Martinotti cell — adapting interneuron targeting layer 1 apical dendrites.
10///
11/// Biophysics: Na+, K+, M-current (Kv7, strong adaptation), T-type Ca2+
12/// (rebound), leak. Overlaps with SST+ phenotype but with stronger
13/// adaptation (higher g_m) and lower rheobase.
14///
15/// Based on Silberberg & Markram 2007 / Toledo-Rodriguez et al. 2005.
16#[derive(Clone, Debug)]
17pub struct MartinottiNeuron {
18    pub v: f64,
19    pub m: f64,
20    pub h: f64,
21    pub n: f64,
22    pub p: f64, // M-current activation
23    pub s: f64, // T-type Ca2+ inactivation
24    // Conductances
25    pub g_na: f64,
26    pub g_k: f64,
27    pub g_m: f64,
28    pub g_t: f64,
29    pub g_l: f64,
30    // Reversal potentials
31    pub e_na: f64,
32    pub e_k: f64,
33    pub e_ca: f64,
34    pub e_l: f64,
35    pub c_m: f64,
36    pub dt: f64,
37    pub v_threshold: f64,
38}
39
40impl MartinottiNeuron {
41    pub fn new() -> Self {
42        Self {
43            v: -65.0,
44            m: 0.02,
45            h: 0.8,
46            n: 0.2,
47            p: 0.0,
48            s: 0.9,
49            g_na: 40.0,
50            g_k: 5.0,
51            g_m: 0.25, // Very strong M-current → pronounced adaptation
52            g_t: 0.01, // T-type Ca2+ (minimal window current)
53            g_l: 0.05, // Leak for resting stability
54            e_na: 50.0,
55            e_k: -90.0,
56            e_ca: 120.0,
57            e_l: -65.0,
58            c_m: 0.8,
59            dt: 0.025,
60            v_threshold: -20.0,
61        }
62    }
63
64    /// Return `[dV, dm, dh, dn, dp, ds]` of the six-state Martinotti system at one
65    /// consistent state. The Na/K activation rates use the L'Hôpital limit at the
66    /// removable Traub-Miles singularity; β_m carries the published `V - V_T - 40`
67    /// offset (an earlier `-17` offset drove the cell into depolarisation block).
68    fn derivatives(
69        &self,
70        v: f64,
71        m: f64,
72        h: f64,
73        n: f64,
74        p: f64,
75        s: f64,
76        current: f64,
77    ) -> [f64; 6] {
78        let dvt = v - (-56.2);
79        let asing = |num: f64, slope: f64, limit: f64| {
80            if num.abs() < 1e-6 {
81                limit
82            } else {
83                num / ((num / slope).exp() - 1.0)
84            }
85        };
86        let alpha_m = -0.32 * asing(dvt - 13.0, -4.0, -4.0);
87        let beta_m = 0.28 * asing(dvt - 40.0, 5.0, 5.0);
88        let alpha_h = 0.128 * (-(dvt - 17.0) / 18.0).exp();
89        let beta_h = 4.0 / (1.0 + (-(dvt - 40.0) / 5.0).exp());
90        let alpha_n = -0.032 * asing(dvt - 15.0, -5.0, -5.0);
91        let beta_n = 0.5 * (-(dvt - 10.0) / 40.0).exp();
92        let dm = alpha_m * (1.0 - m) - beta_m * m;
93        let dh = alpha_h * (1.0 - h) - beta_h * h;
94        let dn = alpha_n * (1.0 - n) - beta_n * n;
95        let p_inf = 1.0 / (1.0 + (-(v + 35.0) / 10.0).exp());
96        let tau_p = 400.0 / (3.3 * ((v + 35.0) / 20.0).exp() + (-(v + 35.0) / 20.0).exp());
97        let dp = (p_inf - p) / tau_p;
98        let m_t_inf = 1.0 / (1.0 + (-(v + 57.0) / 6.2).exp());
99        let s_inf = 1.0 / (1.0 + ((v + 81.0) / 4.0).exp());
100        let tau_s = 30.0 + 200.0 / (1.0 + ((v + 70.0) / 5.0).exp());
101        let ds = (s_inf - s) / tau_s;
102        let i_na = self.g_na * m * m * m * h * (v - self.e_na);
103        let i_k = self.g_k * n * n * n * n * (v - self.e_k);
104        let i_m = self.g_m * p * (v - self.e_k);
105        let i_t = self.g_t * m_t_inf * m_t_inf * s * (v - self.e_ca);
106        let i_l = self.g_l * (v - self.e_l);
107        let dvdt = (-i_na - i_k - i_m - i_t - i_l + current) / self.c_m;
108        [dvdt, dm, dh, dn, dp, ds]
109    }
110
111    /// Return one classical RK4 increment of `[V, m, h, n, p, s]`, holding
112    /// `current` constant across the four stages.
113    fn rk4_substep(&self, st: [f64; 6], current: f64) -> [f64; 6] {
114        let dt = self.dt;
115        let k1 = self.derivatives(st[0], st[1], st[2], st[3], st[4], st[5], current);
116        let mut a = [0.0_f64; 6];
117        for i in 0..6 {
118            a[i] = st[i] + 0.5 * dt * k1[i];
119        }
120        let k2 = self.derivatives(a[0], a[1], a[2], a[3], a[4], a[5], current);
121        for i in 0..6 {
122            a[i] = st[i] + 0.5 * dt * k2[i];
123        }
124        let k3 = self.derivatives(a[0], a[1], a[2], a[3], a[4], a[5], current);
125        for i in 0..6 {
126            a[i] = st[i] + dt * k3[i];
127        }
128        let k4 = self.derivatives(a[0], a[1], a[2], a[3], a[4], a[5], current);
129        let mut out = [0.0_f64; 6];
130        for i in 0..6 {
131            out[i] = st[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
132        }
133        out
134    }
135
136    pub fn step(&mut self, current: f64) -> i32 {
137        let v_prev = self.v;
138        let mut st = [self.v, self.m, self.h, self.n, self.p, self.s];
139        for _ in 0..4 {
140            st = self.rk4_substep(st, current);
141        }
142        self.v = st[0];
143        self.m = st[1];
144        self.h = st[2];
145        self.n = st[3];
146        self.p = st[4];
147        self.s = st[5];
148        if self.v >= self.v_threshold && v_prev < self.v_threshold {
149            1
150        } else {
151            0
152        }
153    }
154
155    pub fn reset(&mut self) {
156        self.v = -65.0;
157        self.m = 0.02;
158        self.h = 0.8;
159        self.n = 0.2;
160        self.p = 0.0;
161        self.s = 0.9;
162    }
163}
164
165impl Default for MartinottiNeuron {
166    fn default() -> Self {
167        Self::new()
168    }
169}
170
171// ═══════════════════════════════════════════════════════════════════
172// Tests
173// ═══════════════════════════════════════════════════════════════════
174
175#[cfg(test)]
176mod tests {
177    use super::super::SSTNeuron;
178    use super::*;
179
180    #[test]
181    fn martinotti_fires_with_input() {
182        let mut n = MartinottiNeuron::new();
183        let spikes: i32 = (0..10000).map(|_| n.step(3.0)).sum();
184        assert!(spikes > 0, "Martinotti must fire with sustained input");
185    }
186
187    #[test]
188    fn martinotti_no_fire_without_input() {
189        let mut n = MartinottiNeuron::new();
190        let spikes: i32 = (0..5000).map(|_| n.step(0.0)).sum();
191        assert_eq!(spikes, 0);
192    }
193
194    #[test]
195    fn martinotti_strong_adaptation() {
196        let mut n = MartinottiNeuron::new();
197        let first: i32 = (0..5000).map(|_| n.step(4.0)).sum();
198        let second: i32 = (0..5000).map(|_| n.step(4.0)).sum();
199        assert!(
200            second <= first + 3,
201            "Martinotti should strongly adapt: first={first}, second={second}"
202        );
203    }
204
205    #[test]
206    fn martinotti_adapts_more_than_sst() {
207        // Martinotti has higher g_m → stronger adaptation
208        let mut mc = MartinottiNeuron::new();
209        let mut sst = SSTNeuron::new();
210        // Use same current magnitude
211        let mc_spikes: i32 = (0..10000).map(|_| mc.step(4.0)).sum();
212        let sst_spikes: i32 = (0..10000).map(|_| sst.step(4.0)).sum();
213        // Martinotti should fire less (stronger adaptation, but lower rheobase too)
214        // At minimum, both should fire
215        assert!(mc_spikes > 0, "Martinotti should fire: got {mc_spikes}");
216        assert!(sst_spikes > 0, "SST should fire: got {sst_spikes}");
217    }
218
219    #[test]
220    fn martinotti_reset_roundtrip() {
221        let mut n = MartinottiNeuron::new();
222        for _ in 0..5000 {
223            n.step(4.0);
224        }
225        n.reset();
226        let mut fresh = MartinottiNeuron::new();
227        let r1: i32 = (0..2000).map(|_| n.step(4.0)).sum();
228        let r2: i32 = (0..2000).map(|_| fresh.step(4.0)).sum();
229        assert_eq!(r1, r2);
230    }
231
232    #[test]
233    fn martinotti_voltage_bounded() {
234        let mut n = MartinottiNeuron::new();
235        for _ in 0..20000 {
236            n.step(10.0);
237        }
238        assert!(n.v.is_finite());
239        assert!(n.p.is_finite());
240    }
241
242    #[test]
243    #[ignore = "wall-clock performance smoke; use Criterion benches for timing evidence"]
244    fn martinotti_performance_10k_steps() {
245        let mut n = MartinottiNeuron::new();
246        let start = std::time::Instant::now();
247        for _ in 0..10_000 {
248            n.step(4.0);
249        }
250        assert!(start.elapsed().as_millis() < 100);
251    }
252}
sc_neurocore_engine/neurons/interneurons/martinotti_neuron.rs

sc_neurocore_engine/neurons/interneurons/
martinotti_neuron.rs
