diff --git a/neurolib/models/wc/loadDefaultParams.py b/neurolib/models/wc/loadDefaultParams.py index 8bf95383..d6d87bff 100644 --- 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 index 5872c0ba..30735884 100644 --- 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 index e2e030b6..974e3ab3 100644 --- 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 index 4d35d6a8..ae213ac1 100644 --- 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)