Differentiale static and time-dependent inputs in WC and WW model (#241)

* differentiate static and time-dependent inputs * update example * revert example change * update aln model with correct input_vars
neurolib-dev · Sep 11, 2023 · 8cdd9ba · 8cdd9ba
1 parent f940967
commit 8cdd9ba
Show file tree

Hide file tree

Showing 4 changed files with 39 additions and 17 deletions.
diff --git a/neurolib/models/wc/loadDefaultParams.py b/neurolib/models/wc/loadDefaultParams.py
@@ -67,8 +67,10 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
     params.mu_inh = 3.0  # inhibitory firing threshold
 
     # values of the external inputs
-    params.exc_ext = 0  # baseline external input to E
-    params.inh_ext = 0  # baseline external input to I
+    params.exc_ext_baseline = 0  # baseline external input to E (static)
+    params.inh_ext_baseline = 0  # baseline external input to I (static)
+    params.exc_ext = 0  # time-dependent external input to E
+    params.inh_ext = 0  # time-dependent external input to I
 
     # ------------------------------------------------------------------------
 

diff --git a/neurolib/models/wc/timeIntegration.py b/neurolib/models/wc/timeIntegration.py
@@ -80,6 +80,9 @@ def timeIntegration(params):
     excs = np.zeros((N, startind + len(t)))
     inhs = np.zeros((N, startind + len(t)))
 
+    exc_ext_baseline = params["exc_ext_baseline"]
+    inh_ext_baseline = params["inh_ext_baseline"]
+
     exc_ext = mu.adjustArrayShape(params["exc_ext"], excs)
     inh_ext = mu.adjustArrayShape(params["inh_ext"], inhs)
 
@@ -125,6 +128,8 @@ def timeIntegration(params):
         inhs,
         exc_input_d,
         inh_input_d,
+        exc_ext_baseline,
+        inh_ext_baseline,
         exc_ext,
         inh_ext,
         tau_exc,
@@ -162,6 +167,8 @@ def timeIntegration_njit_elementwise(
     inhs,
     exc_input_d,
     inh_input_d,
+    exc_ext_baseline,
+    inh_ext_baseline,
     exc_ext,
     inh_ext,
     tau_exc,
@@ -192,10 +199,8 @@ def S_I(x):
         return 1.0 / (1.0 + np.exp(-a_inh * (x - mu_inh)))
 
     for i in range(startind, startind + len(t)):
-
         # loop through all the nodes
         for no in range(N):
-
             # To save memory, noise is saved in the activity array
             noise_exc[no] = excs[no, i]
             noise_inh[no] = inhs[no, i]
@@ -217,8 +222,9 @@ def S_I(x):
                         c_excexc * excs[no, i - 1]  # input from within the excitatory population
                         - c_inhexc * inhs[no, i - 1]  # input from the inhibitory population
                         + exc_input_d[no]  # input from other nodes
-                        + exc_ext[no, i - 1]
-                    )  # external input
+                        + exc_ext_baseline  # baseline external input (static)
+                        + exc_ext[no, i - 1]  # time-dependent external input
+                    )
                     + exc_ou[no]  # ou noise
                 )
             )
@@ -231,8 +237,9 @@ def S_I(x):
                     * S_I(
                         c_excinh * excs[no, i - 1]  # input from the excitatory population
                         - c_inhinh * inhs[no, i - 1]  # input from within the inhibitory population
-                        + inh_ext[no, i - 1]
-                    )  # external input
+                        + inh_ext_baseline  # baseline external input (static)
+                        + inh_ext[no, i - 1]  # time-dependent external input
+                    )
                     + inh_ou[no]  # ou noise
                 )
             )
@@ -246,7 +253,7 @@ def S_I(x):
                 excs[no, i] = 1.0
             if excs[no, i] < 0.0:
                 excs[no, i] = 0.0
-                
+
             if inhs[no, i] > 1.0:
                 inhs[no, i] = 1.0
             if inhs[no, i] < 0.0:

diff --git a/neurolib/models/ww/loadDefaultParams.py b/neurolib/models/ww/loadDefaultParams.py
@@ -65,14 +65,16 @@ def loadDefaultParams(Cmat=None, Dmat=None, seed=None):
     params.tau_exc = 100.0  # ms
     params.gamma_exc = 0.641
     params.w_exc = 1.0
-    params.exc_current = 0.382  # nA
+    params.exc_current_baseline = 0.382  # nA, baseline external input current (static)
+    params.exc_current = 0  # time-dependent external input current to E
 
     params.a_inh = 0.615  # nC^-1
     params.b_inh = 0.177  # kHz
     params.d_inh = 87.0  # ms
     params.tau_inh = 10.0  # ms
     params.w_inh = 0.7
-    params.inh_current = 0.382  # nA
+    params.inh_current_baseline = 0.382  # nA, baseline external input current (static)
+    params.inh_current = 0  # time-dependent external input current to E
 
     params.J_NMDA = 0.15  # nA, excitatory synaptic coupling
     params.J_I = 1.0  # nA, inhibitory synaptic coupling

diff --git a/neurolib/models/ww/timeIntegration.py b/neurolib/models/ww/timeIntegration.py
@@ -25,14 +25,14 @@ def timeIntegration(params):
     tau_exc = params["tau_exc"]
     gamma_exc = params["gamma_exc"]
     w_exc = params["w_exc"]
-    exc_current = params["exc_current"]
+    exc_current_baseline = params["exc_current_baseline"]
 
     a_inh = params["a_inh"]
     b_inh = params["b_inh"]
     d_inh = params["d_inh"]
     tau_inh = params["tau_exc"]
     w_inh = params["w_inh"]
-    inh_current = params["inh_current"]
+    inh_current_baseline = params["inh_current_baseline"]
 
     J_NMDA = params["J_NMDA"]
     J_I = params["J_I"]
@@ -102,6 +102,9 @@ def timeIntegration(params):
     r_exc = np.zeros((N, startind + len(t)))
     r_inh = np.zeros((N, startind + len(t)))
 
+    exc_current = mu.adjustArrayShape(params["exc_current"], r_exc)
+    inh_current = mu.adjustArrayShape(params["inh_current"], r_inh)
+
     # ------------------------------------------------------------------------
     # Set initial values
     # if initial values are just a Nx1 array
@@ -152,12 +155,14 @@ def timeIntegration(params):
         gamma_exc,
         w_exc,
         exc_current,
+        exc_current_baseline,
         a_inh,
         b_inh,
         d_inh,
         tau_inh,
         w_inh,
         inh_current,
+        inh_current_baseline,
         J_NMDA,
         J_I,
         w_ee,
@@ -197,12 +202,14 @@ def timeIntegration_njit_elementwise(
     gamma_exc,
     w_exc,
     exc_current,
+    exc_current_baseline,
     a_inh,
     b_inh,
     d_inh,
     tau_inh,
     w_inh,
     inh_current,
+    inh_current_baseline,
     J_NMDA,
     J_I,
     w_ee,
@@ -232,16 +239,15 @@ def timeIntegration_njit_elementwise(
     r = (a * I - b) / (1.0 - exp(-d * (a * I - b)))
 
     """
+
     #  firing rate transfer function
     def r(I, a, b, d):
         return (a * I - b) / (1.0 - np.exp(-d * (a * I - b)))
 
     ### integrate ODE system:
     for i in range(startind, startind + len(t)):
-
         # loop through all the nodes
         for no in range(N):
-
             # To save memory, noise is saved in the activity array
             noise_se[no] = ses[no, i]
             noise_si[no] = sis[no, i]
@@ -257,8 +263,13 @@ def r(I, a, b, d):
             se = ses[no, i - 1]
             si = sis[no, i - 1]
 
-            I_exc = w_exc * exc_current + w_ee * J_NMDA * se - J_I * si + J_NMDA * ses_input_d[no]
-            I_inh = w_inh * inh_current + J_NMDA * se - si
+            I_exc = (
+                w_exc * (exc_current_baseline + exc_current[no, i - 1])
+                + w_ee * J_NMDA * se
+                - J_I * si
+                + J_NMDA * ses_input_d[no]
+            )
+            I_inh = w_inh * (inh_current_baseline + inh_current[no, i - 1]) + J_NMDA * se - si
 
             r_exc[no, i] = r(I_exc, a_exc, b_exc, d_exc)
             r_inh[no, i] = r(I_inh, a_inh, b_inh, d_inh)