sc_neurocore_engine/neurons/biophysical/
hill_tononi.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 — Hill-Tononi Neuron Model
8
9//! Hill-Tononi thalamocortical sleep and wake dynamics.
10
11/// Hill-Tononi 2005 — thalamocortical sleep/wake with Na-dependent K.
12#[derive(Clone, Debug)]
13pub struct HillTononiNeuron {
14    pub v: f64,
15    pub h_na: f64,
16    pub n_k: f64,
17    pub m_h: f64,
18    pub h_t: f64,
19    pub na_i: f64,
20    pub g_na: f64,
21    pub g_k: f64,
22    pub g_h: f64,
23    pub g_t: f64,
24    pub g_kna: f64,
25    pub g_l: f64,
26    pub e_na: f64,
27    pub e_k: f64,
28    pub e_h: f64,
29    pub e_ca: f64,
30    pub e_l: f64,
31    pub na_pump_max: f64,
32    pub na_eq: f64,
33    pub dt: f64,
34    pub v_threshold: f64,
35}
36
37impl HillTononiNeuron {
38    pub fn new() -> Self {
39        Self {
40            v: -65.0,
41            h_na: 0.6,
42            n_k: 0.3,
43            m_h: 0.0,
44            h_t: 0.9,
45            na_i: 5.0,
46            g_na: 50.0,
47            g_k: 5.0,
48            g_h: 1.0,
49            g_t: 3.0,
50            g_kna: 1.33,
51            g_l: 0.02,
52            e_na: 50.0,
53            e_k: -90.0,
54            e_h: -43.0,
55            e_ca: 120.0,
56            e_l: -70.0,
57            na_pump_max: 20.0,
58            na_eq: 9.5,
59            dt: 0.05,
60            v_threshold: -20.0,
61        }
62    }
63    /// Right-hand side `(dV, dh_na, dn_k, dm_h, dh_t, dna_i)` at one consistent
64    /// state. The sodium and T-type activations are instantaneous; the
65    /// conductance powers use explicit multiplication and the I_KNa Hill exponent
66    /// 3.5 is evaluated as `b*b*b*sqrt(b)` (an IEEE-754 exact decomposition of
67    /// `b.powf(3.5)`) so the Python, Julia, Go, and Mojo backends reproduce the
68    /// trajectory bit-for-bit rather than depending on a per-platform `pow`.
69    fn derivatives(
70        &self,
71        v: f64,
72        h_na: f64,
73        n_k: f64,
74        m_h: f64,
75        h_t: f64,
76        na_i: f64,
77        current: f64,
78    ) -> [f64; 6] {
79        let m_na_inf = 1.0 / (1.0 + (-(v + 38.0) / 10.0).exp());
80        let h_na_inf = 1.0 / (1.0 + ((v + 43.0) / 6.0).exp());
81        let n_k_inf = 1.0 / (1.0 + (-(v + 27.0) / 11.5).exp());
82        let m_h_inf = 1.0 / (1.0 + ((v + 75.0) / 5.5).exp());
83        let m_t_inf = 1.0 / (1.0 + (-(v + 59.0) / 6.2).exp());
84        let h_t_inf = 1.0 / (1.0 + ((v + 83.0) / 4.0).exp());
85        let hill_base = 38.7 / na_i.max(0.01);
86        let w_kna = 0.37 / (1.0 + hill_base * hill_base * hill_base * hill_base.sqrt());
87        let tau_h_na = (1.0 + 10.0 / (1.0 + ((v + 40.0) / 10.0).exp())).max(0.1);
88        let z_n = (v + 50.0) / 25.0;
89        let tau_n_k = (5.0 + 47.0 * (-(z_n * z_n)).exp()).max(0.1);
90        let tau_m_h =
91            (20.0 + 1000.0 / (((v + 71.5) / 14.2).exp() + (-(v + 89.0) / 11.6).exp())).max(1.0);
92        let tau_h_t = if v < -81.0 {
93            (30.8 + 211.4 * ((v + 115.2) / 5.0).exp() / (1.0 + ((v + 86.0) / 3.2).exp())).max(0.1)
94        } else {
95            10.0
96        };
97        let d_h_na = (h_na_inf - h_na) / tau_h_na;
98        let d_n_k = (n_k_inf - n_k) / tau_n_k;
99        let d_m_h = (m_h_inf - m_h) / tau_m_h;
100        let d_h_t = (h_t_inf - h_t) / tau_h_t;
101        let i_na = self.g_na * m_na_inf * m_na_inf * m_na_inf * h_na * (v - self.e_na);
102        let i_k = self.g_k * n_k * n_k * n_k * n_k * (v - self.e_k);
103        let i_h = self.g_h * m_h * (v - self.e_h);
104        let i_t = self.g_t * m_t_inf * m_t_inf * h_t * (v - self.e_ca);
105        let i_kna = self.g_kna * w_kna * (v - self.e_k);
106        let i_l = self.g_l * (v - self.e_l);
107        let d_v = -i_na - i_k - i_h - i_t - i_kna - i_l + current;
108        let d_na_i = -0.001 * i_na - self.na_pump_max * (na_i / (na_i + self.na_eq));
109        [d_v, d_h_na, d_n_k, d_m_h, d_h_t, d_na_i]
110    }
111
112    /// One classical RK4 increment of the six-state vector over `dt`.
113    fn rk4_substep(&self, s: [f64; 6], current: f64) -> [f64; 6] {
114        let dt = self.dt;
115        let k1 = self.derivatives(s[0], s[1], s[2], s[3], s[4], s[5], current);
116        let k2 = self.derivatives(
117            s[0] + 0.5 * dt * k1[0],
118            s[1] + 0.5 * dt * k1[1],
119            s[2] + 0.5 * dt * k1[2],
120            s[3] + 0.5 * dt * k1[3],
121            s[4] + 0.5 * dt * k1[4],
122            s[5] + 0.5 * dt * k1[5],
123            current,
124        );
125        let k3 = self.derivatives(
126            s[0] + 0.5 * dt * k2[0],
127            s[1] + 0.5 * dt * k2[1],
128            s[2] + 0.5 * dt * k2[2],
129            s[3] + 0.5 * dt * k2[3],
130            s[4] + 0.5 * dt * k2[4],
131            s[5] + 0.5 * dt * k2[5],
132            current,
133        );
134        let k4 = self.derivatives(
135            s[0] + dt * k3[0],
136            s[1] + dt * k3[1],
137            s[2] + dt * k3[2],
138            s[3] + dt * k3[3],
139            s[4] + dt * k3[4],
140            s[5] + dt * k3[5],
141            current,
142        );
143        let mut out = [0.0_f64; 6];
144        for i in 0..6 {
145            out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
146        }
147        out
148    }
149
150    pub fn step(&mut self, current: f64) -> i32 {
151        let v_prev = self.v;
152        let s = self.rk4_substep(
153            [self.v, self.h_na, self.n_k, self.m_h, self.h_t, self.na_i],
154            current,
155        );
156        self.v = s[0];
157        self.h_na = s[1];
158        self.n_k = s[2];
159        self.m_h = s[3];
160        self.h_t = s[4];
161        self.na_i = s[5].max(0.0);
162        if self.v >= self.v_threshold && v_prev < self.v_threshold {
163            1
164        } else {
165            0
166        }
167    }
168    pub fn reset(&mut self) {
169        self.v = -65.0;
170        self.h_na = 0.6;
171        self.n_k = 0.3;
172        self.m_h = 0.0;
173        self.h_t = 0.9;
174        self.na_i = 5.0;
175    }
176}
177impl Default for HillTononiNeuron {
178    fn default() -> Self {
179        Self::new()
180    }
181}
182
183#[cfg(test)]
184mod tests {
185    use super::*;
186
187    #[test]
188    fn default_matches_constructor_state() {
189        let default = HillTononiNeuron::default();
190        let constructed = HillTononiNeuron::new();
191        assert_eq!(default.v, constructed.v);
192    }
193
194    #[test]
195    fn hyperpolarized_state_uses_slow_t_current_branch() {
196        let mut n = HillTononiNeuron::new();
197        n.v = -90.0;
198        assert_eq!(n.step(0.0), 0);
199        assert!(n.v.is_finite());
200        assert!(n.h_t.is_finite());
201    }
202
203    #[test]
204    fn hill_tononi_fires() {
205        let mut n = HillTononiNeuron::new();
206        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
207        assert!(t > 0);
208    }
209
210    // -- HillTononi --
211    #[test]
212    fn hill_tononi_silent_without_input() {
213        let mut n = HillTononiNeuron::new();
214        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
215        // May have some spontaneous activity due to Ih
216        assert!(n.v.is_finite());
217    }
218    #[test]
219    fn hill_tononi_reset_clears_state() {
220        let mut n = HillTononiNeuron::new();
221        for _ in 0..100 {
222            n.step(5.0);
223        }
224        n.reset();
225        assert!((n.v - (-65.0)).abs() < 1e-10);
226        assert!((n.na_i - 5.0).abs() < 1e-10);
227    }
228    #[test]
229    fn hill_tononi_extreme_bounded() {
230        let mut n = HillTononiNeuron::new();
231        for _ in 0..200 {
232            n.step(1e4);
233        }
234        assert!(n.v.is_finite());
235    }
236    #[test]
237    fn hill_tononi_na_accumulation() {
238        let mut n = HillTononiNeuron::new();
239        for _ in 0..500 {
240            n.step(5.0);
241        }
242        // Intracellular Na+ should remain finite and non-negative
243        assert!(n.na_i.is_finite());
244        assert!(n.na_i >= 0.0, "Na_i must be non-negative");
245    }
246    #[test]
247    fn hill_tononi_negative_no_crash() {
248        let mut n = HillTononiNeuron::new();
249        for _ in 0..200 {
250            n.step(-10.0);
251        }
252        assert!(n.v.is_finite());
253    }
254    #[test]
255    fn hill_tononi_nan_no_panic() {
256        let mut n = HillTononiNeuron::new();
257        n.step(f64::NAN);
258    }
259}
sc_neurocore_engine/neurons/biophysical/hill_tononi.rs

sc_neurocore_engine/neurons/biophysical/
hill_tononi.rs
