sc_neurocore_engine/neurons/misc/
smooth_muscle.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 — Smooth Muscle Cell Model
8
9//! Calcium-dependent smooth-muscle electrical dynamics.
10
11// ═══════════════════════════════════════════════════════════════════
12// Smooth Muscle Cell
13// ═══════════════════════════════════════════════════════════════════
14
15/// Smooth muscle cell — slow oscillatory electrical activity.
16///
17/// Visceral/vascular smooth muscle with Ca²⁺-dependent oscillations.
18/// Key features distinct from neurons:
19/// - **No fast Na⁺**: depolarisation is Ca²⁺-dependent (L-type)
20/// - **BK (Ca²⁺-activated K⁺)**: repolarisation via BK channels
21/// - **IP3-mediated Ca²⁺ release**: intracellular Ca²⁺ oscillations
22///   from ER/SR via IP3 receptors drive slow waves
23/// - **SERCA pump**: Ca²⁺ reuptake into stores
24/// - **Slow oscillations**: ~3-12 cycles/min (GI slow waves)
25///
26/// dV/dt = (-ICaL - IBK - IL + I_ext) / C_m
27/// dCa/dt = -alpha*ICaL - SERCA(Ca) + IP3_release(Ca, IP3) - Ca/tau_ca
28///
29/// Hirst & Edwards, J Physiol 531:567, 2001.
30/// Imtiaz et al., Biophys J 83:1877, 2002.
31#[derive(Clone, Debug)]
32pub struct SmoothMuscleCell {
33    pub v: f64,
34    pub d: f64,        // CaL activation
35    pub f: f64,        // CaL inactivation
36    pub ca: f64,       // Cytosolic Ca²⁺ (µM)
37    pub ca_store: f64, // ER/SR Ca²⁺ store (µM)
38    pub c_m: f64,
39    pub g_cal: f64, // L-type Ca²⁺
40    pub g_bk: f64,  // BK channel
41    pub g_l: f64,   // Leak
42    pub e_ca: f64,
43    pub e_k: f64,
44    pub e_l: f64,
45    pub tau_ca: f64,   // Ca²⁺ decay (ms)
46    pub v_serca: f64,  // SERCA pump max rate
47    pub k_serca: f64,  // SERCA Km (µM)
48    pub ip3: f64,      // IP3 concentration (µM, constant or input-driven)
49    pub v_ip3r: f64,   // IP3R max release rate
50    pub k_ip3: f64,    // IP3R half-activation (µM)
51    pub k_ca_ip3: f64, // Ca²⁺ co-activation of IP3R (µM)
52    pub kd_bk: f64,    // BK Ca²⁺ half-activation (µM)
53    pub dt: f64,
54    pub sub_steps: usize,
55    pub gain: f64,
56}
57
58impl Default for SmoothMuscleCell {
59    fn default() -> Self {
60        Self::new()
61    }
62}
63
64impl SmoothMuscleCell {
65    pub fn new() -> Self {
66        Self {
67            v: -60.0,
68            d: 0.01,
69            f: 0.95,
70            ca: 0.1,
71            ca_store: 100.0, // ER/SR store (high, ~100 µM)
72            c_m: 1.0,
73            g_cal: 2.0, // L-type Ca²⁺
74            g_bk: 1.0,  // BK
75            g_l: 0.1,
76            e_ca: 60.0,
77            e_k: -80.0,
78            e_l: -50.0,
79            tau_ca: 50.0,
80            v_serca: 0.5,  // SERCA pump rate
81            k_serca: 0.3,  // SERCA Km
82            ip3: 0.5,      // Tonic IP3
83            v_ip3r: 2.0,   // IP3R release rate
84            k_ip3: 0.3,    // IP3R half-act
85            k_ca_ip3: 0.3, // Ca²⁺ co-activation
86            kd_bk: 0.5,    // BK Ca²⁺ Kd
87            dt: 1.0,       // Slow dynamics (1 ms step)
88            sub_steps: 4,
89            gain: 1.0,
90        }
91    }
92
93    #[inline]
94    fn boltz(v: f64, vh: f64, k: f64) -> f64 {
95        1.0 / (1.0 + (-(v - vh) / k).exp())
96    }
97
98    pub fn step(&mut self, current: f64) -> i32 {
99        let input = self.gain * current;
100        let dt_sub = self.dt / self.sub_steps as f64;
101        let v_prev = self.v;
102
103        for _ in 0..self.sub_steps {
104            let v = self.v;
105
106            // CaL d gate (activation)
107            let d_inf = Self::boltz(v, -20.0, 6.0);
108            let tau_d = 5.0 + 20.0 / (1.0 + ((v + 20.0) / 10.0).powi(2)).max(0.01);
109            self.d += dt_sub * (d_inf - self.d) / tau_d;
110
111            // CaL f gate (slow inactivation)
112            let f_inf = Self::boltz(v, -35.0, -8.0);
113            let tau_f = 50.0 + 200.0 / (1.0 + ((v + 35.0) / 10.0).powi(2)).max(0.01);
114            self.f += dt_sub * (f_inf - self.f) / tau_f;
115
116            self.d = self.d.clamp(0.0, 1.0);
117            self.f = self.f.clamp(0.0, 1.0);
118
119            // BK: Ca²⁺-dependent + voltage-dependent
120            let bk_ca = self.ca * self.ca / (self.ca * self.ca + self.kd_bk * self.kd_bk);
121            let bk_v = Self::boltz(v, -10.0, 15.0);
122            let bk_inf = bk_ca * bk_v;
123
124            // Currents
125            let i_cal = self.g_cal * self.d * self.f * (v - self.e_ca);
126            let i_bk = self.g_bk * bk_inf * (v - self.e_k);
127            let i_l = self.g_l * (v - self.e_l);
128
129            let dv = (-(i_cal + i_bk + i_l) + input) / self.c_m;
130            self.v += dt_sub * dv;
131
132            // Ca²⁺ dynamics
133            // Entry via CaL (inward = negative current = Ca²⁺ entry)
134            let ca_entry = if i_cal < 0.0 { -i_cal * 0.01 } else { 0.0 };
135
136            // IP3R release from store: Ca²⁺-induced Ca²⁺ release (CICR)
137            let ip3_act = self.ip3 / (self.ip3 + self.k_ip3);
138            let ca_act = self.ca / (self.ca + self.k_ca_ip3);
139            let ip3_release = self.v_ip3r * ip3_act * ca_act * self.ca_store;
140
141            // SERCA pump (reuptake into store)
142            let serca = self.v_serca * self.ca * self.ca
143                / (self.ca * self.ca + self.k_serca * self.k_serca);
144
145            // Ca²⁺ dynamics
146            self.ca += dt_sub * (ca_entry + ip3_release - serca - self.ca / self.tau_ca);
147            self.ca_store += dt_sub * (serca - ip3_release);
148
149            self.ca = self.ca.max(0.0);
150            self.ca_store = self.ca_store.max(0.0);
151        }
152
153        // Safety
154        self.v = self.v.clamp(-100.0, 40.0);
155        if !self.v.is_finite() {
156            self.v = -60.0;
157        }
158        if !self.ca.is_finite() {
159            self.ca = 0.1;
160        }
161        if !self.ca_store.is_finite() {
162            self.ca_store = 100.0;
163        }
164
165        // "Spike" = slow wave crossing -30 mV
166        if self.v >= -30.0 && v_prev < -30.0 {
167            1
168        } else {
169            0
170        }
171    }
172
173    pub fn reset(&mut self) {
174        *self = Self::new();
175    }
176}
177
178#[cfg(test)]
179mod tests {
180    use super::*;
181
182    // -- Smooth Muscle Cell tests --
183
184    #[test]
185    fn smooth_fires_with_input() {
186        let mut n = SmoothMuscleCell::new();
187        let mut spikes = 0;
188        for _ in 0..10_000 {
189            spikes += n.step(5.0);
190        }
191        assert!(
192            spikes > 0,
193            "Smooth muscle must produce slow waves, got {spikes}"
194        );
195    }
196
197    #[test]
198    fn smooth_silent_without_input() {
199        let mut n = SmoothMuscleCell::new();
200        n.ip3 = 0.0; // No IP3 oscillation driver
201        let mut spikes = 0;
202        for _ in 0..5_000 {
203            spikes += n.step(0.0);
204        }
205        assert!(
206            spikes <= 1,
207            "Should be essentially silent without drive, got {spikes}"
208        );
209    }
210
211    #[test]
212    fn smooth_has_ip3_pathway() {
213        let n = SmoothMuscleCell::new();
214        assert!(n.v_ip3r > 0.0, "Must have IP3R release pathway");
215        assert!(n.ip3 > 0.0, "Must have tonic IP3");
216    }
217
218    #[test]
219    fn smooth_has_serca() {
220        let n = SmoothMuscleCell::new();
221        assert!(n.v_serca > 0.0, "Must have SERCA pump");
222    }
223
224    #[test]
225    fn smooth_ca_store_exists() {
226        let n = SmoothMuscleCell::new();
227        assert!(n.ca_store > 0.0, "Must have ER/SR Ca²⁺ store");
228    }
229
230    #[test]
231    fn smooth_has_bk() {
232        let n = SmoothMuscleCell::new();
233        assert!(n.g_bk > 0.0, "Must have BK (Ca²⁺-activated K) channel");
234    }
235
236    #[test]
237    fn smooth_no_fast_na() {
238        // Smooth muscle does NOT have fast Na channels
239        let n = SmoothMuscleCell::new();
240        // Verify no g_na field exists by checking only CaL and BK
241        assert!(n.g_cal > 0.0, "Depolarisation must be Ca²⁺-dependent");
242    }
243
244    #[test]
245    fn smooth_nan_stays_finite() {
246        let mut n = SmoothMuscleCell::new();
247        n.step(f64::NAN);
248        assert!(n.v.is_finite());
249        assert!(n.ca.is_finite());
250    }
251
252    #[test]
253    fn smooth_reset_clears() {
254        let mut n = SmoothMuscleCell::new();
255        for _ in 0..1000 {
256            n.step(3.0);
257        }
258        n.reset();
259        assert_eq!(n.v, -60.0);
260        assert_eq!(n.ca, 0.1);
261        assert_eq!(n.ca_store, 100.0);
262    }
263
264    #[test]
265    fn smooth_performance_1k_steps() {
266        let start = std::time::Instant::now();
267        let mut n = SmoothMuscleCell::new();
268        for _ in 0..1_000 {
269            std::hint::black_box(n.step(2.0));
270        }
271        let elapsed = start.elapsed();
272        assert!(elapsed.as_millis() < 50, "1k steps must complete in <50ms");
273    }
274
275    #[test]
276    fn smooth_default_matches_constructor() {
277        let default = SmoothMuscleCell::default();
278        let constructed = SmoothMuscleCell::new();
279        assert_eq!(default.v, constructed.v);
280        assert_eq!(default.ca_store, constructed.ca_store);
281        assert_eq!(default.sub_steps, constructed.sub_steps);
282    }
283
284    #[test]
285    fn smooth_nonfinite_calcium_state_recovers_to_baseline() {
286        let mut n = SmoothMuscleCell::new();
287        n.ca = f64::INFINITY;
288        n.ca_store = f64::INFINITY;
289        n.sub_steps = 0;
290        n.step(0.0);
291        assert_eq!(n.ca, 0.1);
292        assert_eq!(n.ca_store, 100.0);
293        assert!(n.v.is_finite());
294    }
295}
sc_neurocore_engine/neurons/misc/smooth_muscle.rs

sc_neurocore_engine/neurons/misc/
smooth_muscle.rs
