From a7f9e79db4b4db79b20c14643cddde01f067a943 Mon Sep 17 00:00:00 2001
From: Alec Hoyland <alec.hoyland@posteo.net>
Date: Mon, 13 May 2019 15:27:35 -0400
Subject: [PATCH] Spin off `u` into a separate variable

---
 c++/synapses/borgers/NMDAergic.hpp | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/c++/synapses/borgers/NMDAergic.hpp b/c++/synapses/borgers/NMDAergic.hpp
index 0ce949a8..4731c3cb 100644
--- a/c++/synapses/borgers/NMDAergic.hpp
+++ b/c++/synapses/borgers/NMDAergic.hpp
@@ -16,21 +16,24 @@ class NMDAergic: public synapse {
 public:
 
   double Mg; // Mg++ concentration in millimolar (mM)
-  double tau_r;
-  double tau_d;
+  double u; // Mg++ block gating variable
+  double tau_r; // rising-phase time constant (ms)
+  double tau_d; // decaying-phase time constant (ms)
 
     // specify parameters + initial conditions
-    NMDAergic(double g_, double s_, double E_, double Mg_, double tau_r_, double tau_d_)
+    NMDAergic(double g_, double s_, double u_, double E_, double Mg_, double tau_r_, double tau_d_)
     {
         gmax = g_;
         E = E_;
         s = s_;
+        u = u_;
         Mg = Mg_;
         tau_r = tau_r_;
         tau_d = tau_d_;
 
         // defaults
         if (isnan (s)) { s = 0; }
+        if (isnan(u)) { u = 0; }
         if (isnan (gmax)) { gmax = 0; }
         if (isnan (E)) { E = 0; }
         if (isnan (Mg)) { Mg = 1; }
@@ -60,7 +63,7 @@ class NMDAergic: public synapse {
 
 int NMDAergic::getFullStateSize()
 {
-    return 2;
+    return 3;
 }
 
 // // //
@@ -93,8 +96,9 @@ void NMDAergic::integrate(void)
     double sinf = s_inf(ss);
 
     // integrate using exponential Euler
-    s = (sinf + (s - sinf)*exp(-dt/tau_s(ss))) * 1.0 / (1.0 + Mg / 3.57 * exp(-0.062*V_post));
-    g = gmax*s;
+    s = (sinf + (s - sinf)*exp(-dt/tau_s(ss)));
+    u = 1.0 / (1.0 + Mg / 3.57 * exp(-0.062*V_post));
+    g = gmax*s*u;
 }
 
 void NMDAergic::integrateMS(int k, double V, double Ca)
@@ -132,7 +136,7 @@ void NMDAergic::integrateMS(int k, double V, double Ca)
         // last step
 
         s = s + (k_s[0] + 2*k_s[1] + 2*k_s[2] + k_s[3])/6;
-        s = s * 1.0 / (1.0 + Mg / 3.57 * exp(-0.062*V_post));
+        u = 1.0 / (1.0 + Mg / 3.57 * exp(-0.062*V_post));
 
         if (s < 0) {s = 0;}
         if (s > 1) {s = 1;}
@@ -145,11 +149,14 @@ int NMDAergic::getFullState(double *syn_state, int idx)
 {
     // give it the current synapse variable
     syn_state[idx] = s;
-
+    idx++;
+    
+    // give it the current magnesium block variable
+    syn_state[idx] = u;
     idx++;
 
     // also return the current from this synapse
-    syn_state[idx] = gmax*s*(post_syn->V - E);
+    syn_state[idx] = gmax*s*u*(post_syn->V - E);
     idx++;
     return idx;
 }