sc_neurocore_engine/neurons/misc/
cardiac_purkinje.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 — Cardiac Purkinje Fibre Model
8
9//! DiFrancesco-Noble cardiac Purkinje fibre dynamics.
10
11// ═══════════════════════════════════════════════════════════════════
12// Cardiac Purkinje Fibre
13// ═══════════════════════════════════════════════════════════════════
14
15/// Cardiac Purkinje fibre — DiFrancesco-Noble 1985 model.
16///
17/// Specialised cardiac conduction cell with long action potentials
18/// (~300 ms) and pacemaker capability via If (funny/HCN current).
19///
20/// 6 major ionic currents:
21/// - **INa** (fast Na, m³h): rapid depolarisation (phase 0)
22/// - **ICaL** (L-type Ca²⁺, d·f): plateau maintenance (phase 2)
23/// - **IKr** (rapid delayed rectifier K, x_r): phase 3 repolarisation
24/// - **IK1** (inward rectifier K): resting potential stabilisation
25/// - **If** (funny current, HCN, y): pacemaker depolarisation (phase 4)
26/// - **IL** (leak)
27///
28/// Action potential phases:
29/// 0 — rapid depolarisation (INa)
30/// 1 — early repolarisation notch
31/// 2 — plateau (ICaL vs IKr balance)
32/// 3 — repolarisation (IKr dominates)
33/// 4 — pacemaker depolarisation (If)
34///
35/// Uses 10 sub-steps (dt_sub = 0.05 ms) for Na gating stability.
36///
37/// DiFrancesco & Noble, Phil Trans R Soc Lond B 307:353, 1985.
38/// Noble, J Physiol 353:1, 1984 (review).
39#[derive(Clone, Debug)]
40pub struct CardiacPurkinjeFibre {
41    pub v: f64,
42    pub m: f64,   // Na activation
43    pub h: f64,   // Na inactivation
44    pub d: f64,   // CaL activation
45    pub f: f64,   // CaL inactivation
46    pub x_r: f64, // IKr activation
47    pub y: f64,   // If (HCN) activation
48    pub c_m: f64,
49    pub g_na: f64,
50    pub g_cal: f64,
51    pub g_kr: f64,
52    pub g_k1: f64,
53    pub g_f: f64, // Funny current conductance
54    pub g_l: f64,
55    pub e_na: f64,
56    pub e_ca: f64,
57    pub e_k: f64,
58    pub e_f: f64, // If reversal (~-20 mV, mixed cation)
59    pub e_l: f64,
60    pub dt: f64,
61    pub sub_steps: usize,
62    pub gain: f64,
63}
64
65impl Default for CardiacPurkinjeFibre {
66    fn default() -> Self {
67        Self::new()
68    }
69}
70
71impl CardiacPurkinjeFibre {
72    pub fn new() -> Self {
73        Self {
74            v: -85.0,
75            m: 0.001,
76            h: 0.99,
77            d: 0.001,
78            f: 0.99,
79            x_r: 0.01,
80            y: 0.05,
81            c_m: 1.0,
82            g_na: 15.0,  // Fast Na
83            g_cal: 0.05, // L-type Ca²⁺ (small but sustains plateau)
84            g_kr: 0.015, // Rapid delayed rectifier
85            g_k1: 0.4,   // Inward rectifier
86            g_f: 0.01,   // Funny current (pacemaker)
87            g_l: 0.03,
88            e_na: 40.0,
89            e_ca: 65.0,
90            e_k: -90.0,
91            e_f: -20.0, // Mixed Na⁺/K⁺ cation
92            e_l: -50.0,
93            dt: 0.5,
94            sub_steps: 10, // dt_sub = 0.05 ms
95            gain: 1.0,
96        }
97    }
98
99    #[inline]
100    fn boltz(v: f64, vh: f64, k: f64) -> f64 {
101        1.0 / (1.0 + (-(v - vh) / k).exp())
102    }
103
104    pub fn step(&mut self, current: f64) -> i32 {
105        let input = self.gain * current;
106        let dt_sub = self.dt / self.sub_steps as f64;
107        let v_prev = self.v;
108
109        for _ in 0..self.sub_steps {
110            let v = self.v;
111
112            // Na m gate (fast)
113            let m_inf = Self::boltz(v, -40.0, 8.0);
114            let tau_m = 0.05 + 0.3 / (1.0 + ((v + 40.0) / 10.0).powi(2)).max(0.01);
115            self.m += dt_sub * (m_inf - self.m) / tau_m;
116
117            // Na h gate (inactivation)
118            let h_inf = Self::boltz(v, -65.0, -7.0);
119            let tau_h = 0.5 + 8.0 / (1.0 + ((v + 65.0) / 15.0).powi(2)).max(0.01);
120            self.h += dt_sub * (h_inf - self.h) / tau_h;
121
122            // CaL d gate (activation)
123            let d_inf = Self::boltz(v, -10.0, 6.0);
124            let tau_d = 2.0 + 5.0 / (1.0 + ((v + 10.0) / 10.0).powi(2)).max(0.01);
125            self.d += dt_sub * (d_inf - self.d) / tau_d;
126
127            // CaL f gate (inactivation, slow)
128            let f_inf = Self::boltz(v, -30.0, -8.0);
129            let tau_f = 20.0 + 100.0 / (1.0 + ((v + 30.0) / 10.0).powi(2)).max(0.01);
130            self.f += dt_sub * (f_inf - self.f) / tau_f;
131
132            // IKr x_r gate (slow activation)
133            let xr_inf = Self::boltz(v, -20.0, 10.0);
134            let tau_xr = 50.0 + 200.0 / (1.0 + ((v + 20.0) / 15.0).powi(2)).max(0.01);
135            self.x_r += dt_sub * (xr_inf - self.x_r) / tau_xr;
136
137            // If y gate (activates at hyperpolarised V)
138            let y_inf = Self::boltz(v, -80.0, -10.0);
139            let tau_y = 100.0 + 500.0 / (1.0 + ((v + 80.0) / 20.0).powi(2)).max(0.01);
140            self.y += dt_sub * (y_inf - self.y) / tau_y;
141
142            // Clamp gates
143            self.m = self.m.clamp(0.0, 1.0);
144            self.h = self.h.clamp(0.0, 1.0);
145            self.d = self.d.clamp(0.0, 1.0);
146            self.f = self.f.clamp(0.0, 1.0);
147            self.x_r = self.x_r.clamp(0.0, 1.0);
148            self.y = self.y.clamp(0.0, 1.0);
149
150            // IK1: inward rectifier (voltage-dependent, Boltzmann)
151            let k1_inf = 1.0 / (1.0 + ((v - self.e_k + 10.0) / 10.0).exp());
152
153            // Currents
154            let i_na = self.g_na * self.m.powi(3) * self.h * (v - self.e_na);
155            let i_cal = self.g_cal * self.d * self.f * (v - self.e_ca);
156            let i_kr = self.g_kr * self.x_r * (v - self.e_k);
157            let i_k1 = self.g_k1 * k1_inf * (v - self.e_k);
158            let i_f = self.g_f * self.y * (v - self.e_f);
159            let i_l = self.g_l * (v - self.e_l);
160
161            let dv = (-(i_na + i_cal + i_kr + i_k1 + i_f + i_l) + input) / self.c_m;
162            self.v += dt_sub * dv;
163        }
164
165        // Safety
166        self.v = self.v.clamp(-120.0, 60.0);
167        if !self.v.is_finite() {
168            self.v = -85.0;
169        }
170        if !self.m.is_finite() {
171            self.m = 0.001;
172        }
173        if !self.h.is_finite() {
174            self.h = 0.99;
175        }
176        if !self.d.is_finite() {
177            self.d = 0.001;
178        }
179        if !self.f.is_finite() {
180            self.f = 0.99;
181        }
182        if !self.x_r.is_finite() {
183            self.x_r = 0.01;
184        }
185        if !self.y.is_finite() {
186            self.y = 0.05;
187        }
188
189        // Spike: V crosses -20 mV upward
190        if self.v >= -20.0 && v_prev < -20.0 {
191            1
192        } else {
193            0
194        }
195    }
196
197    pub fn reset(&mut self) {
198        *self = Self::new();
199    }
200}
201
202#[cfg(test)]
203mod tests {
204    use super::*;
205
206    // -- Cardiac Purkinje Fibre tests --
207
208    #[test]
209    fn cardiac_fires_with_input() {
210        let mut n = CardiacPurkinjeFibre::new();
211        let mut spikes = 0;
212        for _ in 0..2_000 {
213            spikes += n.step(5.0);
214        }
215        assert!(
216            spikes > 0,
217            "Cardiac Purkinje must fire with input, got {spikes}"
218        );
219    }
220
221    #[test]
222    fn cardiac_silent_without_input() {
223        // Without If-driven pacemaking (test with g_f=0)
224        let mut n = CardiacPurkinjeFibre::new();
225        n.g_f = 0.0; // Disable pacemaker
226        let mut spikes = 0;
227        for _ in 0..5_000 {
228            spikes += n.step(0.0);
229        }
230        assert!(
231            spikes <= 1,
232            "Must be essentially silent without pacemaker, got {spikes}"
233        );
234    }
235
236    #[test]
237    fn cardiac_has_funny_current() {
238        // If (HCN) is the hallmark pacemaker current
239        let n = CardiacPurkinjeFibre::new();
240        assert!(n.g_f > 0.0, "Must have funny current (If/HCN)");
241    }
242
243    #[test]
244    fn cardiac_has_cal() {
245        // L-type Ca²⁺ sustains the plateau
246        let n = CardiacPurkinjeFibre::new();
247        assert!(n.g_cal > 0.0, "Must have L-type Ca²⁺ for plateau");
248    }
249
250    #[test]
251    fn cardiac_has_inward_rectifier() {
252        // IK1 stabilises resting potential
253        let n = CardiacPurkinjeFibre::new();
254        assert!(n.g_k1 > 0.0, "Must have IK1 inward rectifier");
255    }
256
257    #[test]
258    fn cardiac_six_currents() {
259        let n = CardiacPurkinjeFibre::new();
260        assert!(
261            n.g_na > 0.0
262                && n.g_cal > 0.0
263                && n.g_kr > 0.0
264                && n.g_k1 > 0.0
265                && n.g_f > 0.0
266                && n.g_l > 0.0,
267            "Must have all 6 currents"
268        );
269    }
270
271    #[test]
272    fn cardiac_gating_evolves() {
273        let mut n = CardiacPurkinjeFibre::new();
274        let d0 = n.d;
275        let y0 = n.y;
276        for _ in 0..200 {
277            n.step(5.0);
278        }
279        assert!(n.d != d0 || n.y != y0, "Gating must evolve");
280    }
281
282    #[test]
283    fn cardiac_nan_input_stays_finite() {
284        let mut n = CardiacPurkinjeFibre::new();
285        n.step(f64::NAN);
286        assert!(n.v.is_finite());
287    }
288
289    #[test]
290    fn cardiac_reset_clears_state() {
291        let mut n = CardiacPurkinjeFibre::new();
292        for _ in 0..500 {
293            n.step(5.0);
294        }
295        n.reset();
296        assert_eq!(n.v, -85.0);
297        assert_eq!(n.m, 0.001);
298    }
299
300    #[test]
301    fn cardiac_performance_1k_steps() {
302        let start = std::time::Instant::now();
303        let mut n = CardiacPurkinjeFibre::new();
304        for _ in 0..1_000 {
305            std::hint::black_box(n.step(3.0));
306        }
307        let elapsed = start.elapsed();
308        assert!(elapsed.as_millis() < 50, "1k steps must complete in <50ms");
309    }
310
311    #[test]
312    fn cardiac_default_matches_constructor() {
313        let default = CardiacPurkinjeFibre::default();
314        let constructed = CardiacPurkinjeFibre::new();
315        assert_eq!(default.v, constructed.v);
316        assert_eq!(default.g_f, constructed.g_f);
317        assert_eq!(default.sub_steps, constructed.sub_steps);
318    }
319}
sc_neurocore_engine/neurons/misc/cardiac_purkinje.rs

sc_neurocore_engine/neurons/misc/
cardiac_purkinje.rs
