sc_neurocore_engine/neurons/biophysical/
durstewitz_dopamine.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 — Durstewitz Dopamine Neuron Model
8
9//! Durstewitz prefrontal-cortex neuron dynamics with D1 modulation.
10
11/// Durstewitz PFC neuron with D1 dopamine modulation. Durstewitz et al. 2000.
12#[derive(Clone, Debug)]
13pub struct DurstewitzDopamineNeuron {
14    pub v: f64,
15    pub h_na: f64,
16    pub n_k: f64,
17    pub g_na: f64,
18    pub g_k: f64,
19    pub g_nmda: f64,
20    pub g_l: f64,
21    pub e_na: f64,
22    pub e_k: f64,
23    pub e_nmda: f64,
24    pub e_l: f64,
25    pub mg: f64,
26    pub d1_level: f64,
27    pub g_nmda_scale: f64,
28    pub g_k_scale: f64,
29    pub v_shift_na: f64,
30    pub dt: f64,
31    pub v_threshold: f64,
32}
33
34impl DurstewitzDopamineNeuron {
35    pub fn new() -> Self {
36        Self {
37            v: -65.0,
38            h_na: 0.7,
39            n_k: 0.2,
40            g_na: 45.0,
41            g_k: 18.0,
42            g_nmda: 0.5,
43            g_l: 0.02,
44            e_na: 55.0,
45            e_k: -80.0,
46            e_nmda: 0.0,
47            e_l: -65.0,
48            mg: 1.0,
49            d1_level: 0.0,
50            g_nmda_scale: 2.5,
51            g_k_scale: 1.5,
52            v_shift_na: -5.0,
53            dt: 0.05,
54            v_threshold: -20.0,
55        }
56    }
57    /// Right-hand side ``(dV, dh_na, dn_k)`` evaluated from one consistent state.
58    ///
59    /// The sodium activation ``m_∞`` is instantaneous, so it is recomputed from
60    /// `v` at every RK4 stage. The conductance powers use explicit multiplication
61    /// and the Mg²⁺ block keeps the `mg / 3.57 * exp` operand order so the
62    /// Python, Julia, Go, and Mojo backends reproduce the trajectory bit-for-bit.
63    fn derivatives(&self, v: f64, h_na: f64, n_k: f64, current: f64) -> [f64; 3] {
64        let v_sh = self.d1_level * self.v_shift_na;
65        let m_na_inf = 1.0 / (1.0 + (-(v + 30.0 + v_sh) / 9.5).exp());
66        let h_na_inf = 1.0 / (1.0 + ((v + 53.0) / 7.0).exp());
67        let n_k_inf = 1.0 / (1.0 + (-(v + 30.0) / 10.0).exp());
68        let tau_h = 0.5 + 14.0 / (1.0 + ((v + 50.0) / 12.0).exp());
69        let tau_n = 1.0 + 11.0 / (1.0 + ((v + 40.0) / 10.0).exp());
70        let d_h_na = (h_na_inf - h_na) / tau_h;
71        let d_n_k = (n_k_inf - n_k) / tau_n;
72        let mg_block = 1.0 / (1.0 + self.mg / 3.57 * (-0.062 * v).exp());
73        let nmda_g = self.g_nmda * (1.0 + self.d1_level * (self.g_nmda_scale - 1.0));
74        let k_g = self.g_k * (1.0 + self.d1_level * (self.g_k_scale - 1.0));
75        let i_na = self.g_na * m_na_inf * m_na_inf * m_na_inf * h_na * (v - self.e_na);
76        let i_k = k_g * n_k * n_k * n_k * n_k * (v - self.e_k);
77        let i_nmda = nmda_g * mg_block * (v - self.e_nmda);
78        let i_l = self.g_l * (v - self.e_l);
79        let d_v = -i_na - i_k - i_nmda - i_l + current;
80        [d_v, d_h_na, d_n_k]
81    }
82
83    /// One classical RK4 increment of the `(V, h_na, n_k)` vector over `dt`.
84    fn rk4_substep(&self, s: [f64; 3], current: f64) -> [f64; 3] {
85        let dt = self.dt;
86        let k1 = self.derivatives(s[0], s[1], s[2], current);
87        let k2 = self.derivatives(
88            s[0] + 0.5 * dt * k1[0],
89            s[1] + 0.5 * dt * k1[1],
90            s[2] + 0.5 * dt * k1[2],
91            current,
92        );
93        let k3 = self.derivatives(
94            s[0] + 0.5 * dt * k2[0],
95            s[1] + 0.5 * dt * k2[1],
96            s[2] + 0.5 * dt * k2[2],
97            current,
98        );
99        let k4 = self.derivatives(
100            s[0] + dt * k3[0],
101            s[1] + dt * k3[1],
102            s[2] + dt * k3[2],
103            current,
104        );
105        let mut out = [0.0_f64; 3];
106        for i in 0..3 {
107            out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
108        }
109        out
110    }
111
112    pub fn step(&mut self, current: f64) -> i32 {
113        let v_prev = self.v;
114        let s = self.rk4_substep([self.v, self.h_na, self.n_k], current);
115        self.v = s[0];
116        self.h_na = s[1];
117        self.n_k = s[2];
118        if self.v >= self.v_threshold && v_prev < self.v_threshold {
119            1
120        } else {
121            0
122        }
123    }
124    pub fn reset(&mut self) {
125        self.v = -65.0;
126        self.h_na = 0.7;
127        self.n_k = 0.2;
128    }
129}
130impl Default for DurstewitzDopamineNeuron {
131    fn default() -> Self {
132        Self::new()
133    }
134}
135
136#[cfg(test)]
137mod tests {
138    use super::*;
139
140    #[test]
141    fn default_matches_constructor_state() {
142        let default = DurstewitzDopamineNeuron::default();
143        let constructed = DurstewitzDopamineNeuron::new();
144        assert_eq!(default.v, constructed.v);
145    }
146
147    #[test]
148    fn durstewitz_fires() {
149        let mut n = DurstewitzDopamineNeuron::new();
150        let t: i32 = (0..1000).map(|_| n.step(3.0)).sum();
151        assert!(t > 0);
152    }
153
154    // -- DurstewitzDopamine --
155    #[test]
156    fn durstewitz_low_activity_zero_input() {
157        let mut n = DurstewitzDopamineNeuron::new();
158        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
159        // NMDA tonic conductance can produce spontaneous activity
160        assert!(n.v.is_finite());
161    }
162    #[test]
163    fn durstewitz_reset_clears_state() {
164        let mut n = DurstewitzDopamineNeuron::new();
165        for _ in 0..100 {
166            n.step(3.0);
167        }
168        n.reset();
169        assert!((n.v - (-65.0)).abs() < 1e-10);
170    }
171    #[test]
172    fn durstewitz_extreme_bounded() {
173        let mut n = DurstewitzDopamineNeuron::new();
174        for _ in 0..200 {
175            n.step(1e4);
176        }
177        assert!(n.v.is_finite());
178    }
179    #[test]
180    fn durstewitz_d1_modulation() {
181        // D1 dopamine should increase NMDA and shift Na activation
182        let mut n_d1 = DurstewitzDopamineNeuron::new();
183        n_d1.d1_level = 1.0;
184        let mut n_no = DurstewitzDopamineNeuron::new();
185        n_no.d1_level = 0.0;
186        for _ in 0..1000 {
187            n_d1.step(3.0);
188        }
189        for _ in 0..1000 {
190            n_no.step(3.0);
191        }
192        // Both should remain stable; D1 changes effective conductances
193        assert!(n_d1.v.is_finite() && n_no.v.is_finite());
194    }
195    #[test]
196    fn durstewitz_mg_block() {
197        let n = DurstewitzDopamineNeuron::new();
198        // At rest (-65 mV), Mg²⁺ block should be high
199        let block = 1.0 / (1.0 + n.mg * (-0.062 * n.v).exp() / 3.57);
200        assert!(block < 0.1, "Mg²⁺ block at rest should be high: {}", block);
201    }
202    #[test]
203    fn durstewitz_negative_no_crash() {
204        let mut n = DurstewitzDopamineNeuron::new();
205        for _ in 0..200 {
206            n.step(-10.0);
207        }
208        assert!(n.v.is_finite());
209    }
210}
sc_neurocore_engine/neurons/biophysical/durstewitz_dopamine.rs

sc_neurocore_engine/neurons/biophysical/
durstewitz_dopamine.rs
