sc_neurocore_engine/neurons/simple_spiking/
chay.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 — Chay Neuron Model
8
9//! Chay pancreatic beta-cell dynamics.
10
11/// Chay 1985 — pancreatic beta cell with Ca-dependent K.
12#[derive(Clone, Debug)]
13pub struct ChayNeuron {
14    pub v: f64,
15    pub n: f64,
16    pub ca: f64,
17    pub g_ca: f64,
18    pub g_k: f64,
19    pub g_kca: f64,
20    pub g_l: f64,
21    pub e_ca: f64,
22    pub e_k: f64,
23    pub e_l: f64,
24    pub rho: f64,
25    pub alpha_ca: f64,
26    pub k_ca: f64,
27    pub dt: f64,
28    pub v_threshold: f64,
29}
30
31impl ChayNeuron {
32    pub fn new() -> Self {
33        Self {
34            v: -50.0,
35            n: 0.1,
36            ca: 0.1,
37            g_ca: 25.0,
38            g_k: 1400.0,
39            g_kca: 12.0,
40            g_l: 7.0,
41            e_ca: 100.0,
42            e_k: -75.0,
43            e_l: -40.0,
44            rho: 0.00015,
45            alpha_ca: 0.002,
46            k_ca: 0.04,
47            dt: 0.02,
48            v_threshold: -20.0,
49        }
50    }
51    pub fn step(&mut self, current: f64) -> i32 {
52        if !current.is_finite() || !self.dt.is_finite() || self.dt <= 0.0 {
53            return 0;
54        }
55
56        let v_initial = self.v;
57        let mut v = self.v;
58        let mut n = self.n;
59        let mut ca = self.ca;
60        let substeps = (self.dt / 0.001_f64).ceil().max(1.0) as usize;
61        let h = self.dt / substeps as f64;
62        let mut crossed = false;
63
64        for _ in 0..substeps {
65            let m_inf = 1.0 / (1.0 + (-(v + 25.0) / 8.0).clamp(-700.0, 700.0).exp());
66            let n_inf = 1.0 / (1.0 + (-(v + 18.0) / 14.0).clamp(-700.0, 700.0).exp());
67            let d = (v + 18.0).abs().max(0.01);
68            let tau_n = 1.0 / (0.01 * d);
69            let ca_denominator = ca + 1.0;
70            if ca_denominator <= 0.0 {
71                return 0;
72            }
73            let kca_act = ca / ca_denominator;
74            let i_ca = self.g_ca * m_inf * (v - self.e_ca);
75            let i_k = self.g_k * n * (v - self.e_k);
76            let i_kca = self.g_kca * kca_act * (v - self.e_k);
77            let i_l = self.g_l * (v - self.e_l);
78
79            let v_next = v + (-i_ca - i_k - i_kca - i_l + current) * h;
80            let n_next = n + (n_inf - n) / tau_n.max(0.01) * h;
81            let ca_next = ca + self.rho * (-self.alpha_ca * i_ca - self.k_ca * ca) * h;
82            if !v_next.is_finite()
83                || !n_next.is_finite()
84                || !ca_next.is_finite()
85                || !(-200.0..=200.0).contains(&v_next)
86                || !(0.0..=1.0).contains(&n_next)
87                || !(0.0..=100.0).contains(&ca_next)
88            {
89                return 0;
90            }
91            crossed = crossed || (v_next >= self.v_threshold && v < self.v_threshold);
92            v = v_next;
93            n = n_next;
94            ca = ca_next;
95        }
96
97        self.v = v;
98        self.n = n;
99        self.ca = ca;
100        if crossed && v_initial < self.v_threshold {
101            1
102        } else {
103            0
104        }
105    }
106    pub fn reset(&mut self) {
107        self.v = -50.0;
108        self.n = 0.1;
109        self.ca = 0.1;
110    }
111}
112impl Default for ChayNeuron {
113    fn default() -> Self {
114        Self::new()
115    }
116}
117
118#[cfg(test)]
119mod tests {
120    use super::*;
121
122    #[test]
123    fn default_matches_constructor_state() {
124        let default = ChayNeuron::default();
125        let constructed = ChayNeuron::new();
126        assert_eq!(default.v, constructed.v);
127    }
128
129    #[test]
130    fn chay_drive_changes_state_without_leaving_physical_bounds() {
131        let mut rest = ChayNeuron::new();
132        let mut driven = ChayNeuron::new();
133        for _ in 0..500 {
134            rest.step(0.0);
135            driven.step(5.0);
136        }
137        assert!(driven.v > rest.v);
138        assert!((0.0..=1.0).contains(&driven.n));
139        assert!(driven.ca >= 0.0);
140    }
141
142    #[test]
143    fn chay_reset_clears_state() {
144        let mut n = ChayNeuron::new();
145        for _ in 0..1000 {
146            n.step(20.0);
147        }
148        n.reset();
149        assert!((n.v - (-50.0)).abs() < 1e-10);
150    }
151
152    #[test]
153    fn chay_bounded() {
154        let mut n = ChayNeuron::new();
155        for _ in 0..5000 {
156            n.step(200.0);
157        }
158        assert!(n.v.is_finite());
159    }
160
161    #[test]
162    fn chay_ca_nonneg() {
163        let mut n = ChayNeuron::new();
164        for _ in 0..5000 {
165            n.step(20.0);
166        }
167        assert!(n.ca >= 0.0, "Ca²⁺ must be non-negative");
168    }
169
170    #[test]
171    fn chay_nan_no_panic() {
172        ChayNeuron::new().step(f64::NAN);
173    }
174
175    #[test]
176    fn chay_negative_no_crash() {
177        let mut n = ChayNeuron::new();
178        for _ in 0..500 {
179            n.step(-10.0);
180        }
181        assert!(n.v.is_finite());
182    }
183}
sc_neurocore_engine/neurons/simple_spiking/chay.rs

sc_neurocore_engine/neurons/simple_spiking/
chay.rs
