Created
May 15, 2025 09:35
-
-
Save mstimberg/4195eab8d71e99c28b58a089cea58392 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from brian2 import * | |
| # parameter dictionary definition | |
| params = { | |
| 'defaultclock.dt' : '0.1*ms', | |
| 'Init_time' : '50000*ms', | |
| 'buffer_time' : '1200*ms', | |
| 'run_time' : '1200*ms', | |
| 'Cm' : '1* uF/cm**2', | |
| 'S' : '20e3 * um**2', | |
| 'gL' : '.05 * mS/cm**2', | |
| 'EL' : '-60 * mV', | |
| 'Vr' : '-60 * mV', | |
| 'Vt' : '-50 * mV', | |
| 'DeltaT' : '2.5 * mV', | |
| 'tauw' : '600 * ms', | |
| 'N_i' : '150', | |
| 'N_e' : '1600', | |
| 'N_t' : '200', | |
| 'N_r' : '100', | |
| 'g_E_T' : '0.15 *mS/ cm**2', | |
| 'g_E_E' : '0.2 *mS/ cm**2', | |
| 'g_E_I' : '0.7 *mS/ cm**2', | |
| 'g_I_T' : '0.2 *mS/ cm**2', | |
| 'g_I_E' : '0.6 *mS/ cm**2', | |
| 'g_I_I' : '0.55 *mS/ cm**2', | |
| 'g_R_E' : '0.35 *mS/ cm**2', | |
| 'g_T_E' : '0.35 *mS/ cm**2', | |
| 'g_R_T' : '0.085 *mS/ cm**2', | |
| 'g_T_R' : '2.84 *mS/ cm**2', | |
| 'g_R_R' : '2.84 *mS/ cm**2', | |
| 'E_Ampa' : '0.0*mV', | |
| 'E_Gaba' : '-80.0*mV', | |
| 'tau_Ampa' : '2.0 *ms', | |
| 'tau_Gaba' : '3.0 *ms', | |
| 'Kappa_e_t' : '100.0', | |
| 'Kappa_i_t' : '75.0', | |
| 'Kappa_e_e' : '200.0', | |
| 'Kappa_i_e' : '400.0', | |
| 'Kappa_e_i' : '25.0', | |
| 'Kappa_i_i' : '25.0', | |
| 'Kappa_t_e' : '32.0', | |
| 'Kappa_r_e' : '32.0', | |
| 'Kappa_r_t' : '2.0', | |
| 'Kappa_t_r' : '8.0', | |
| 'Kappa_r_r' : '8.0', | |
| 'syn_weight' : '1', | |
| 'Lambda' : '0.1 * mm', | |
| 'R_Cx' : '0.20 * mm', | |
| 'R_Th' : '0.10 * mm', | |
| 'Po_N_exc' : '100', | |
| 'Po_N_inh' : '25', | |
| 'Po_Rate': '25*Hz', | |
| 'Po_Amp_T': '0.2*mV', | |
| 'Po_Amp_E': '0.2*mV', | |
| 'Po_Amp_I': '0.07*mV', | |
| 'Po_Amp_R': '0.2*mV', | |
| } | |
| def setParams(d): | |
| '''Sets a dictionary of parameters as global python variables. | |
| Args: | |
| d (dict): The dictionary of parameters (param_name, param_vals) | |
| ''' | |
| for param, val in d.items(): | |
| exec(param + '=' + val, globals()) | |
| def setParam(param, val, d): | |
| '''Sets a single parameter/value pair as a global python variable and adds it to data export dict. | |
| Args: | |
| param (str): The name of the parameter variable | |
| val (str): The value of the parameter variable | |
| d (dict): The dictionary of parameters (param_name, param_vals) | |
| ''' | |
| exec(param + '=' + val, globals()) | |
| d[param] = val | |
| setParams(params) | |
| #AdEx IF Model | |
| eqs_aeIF = """ | |
| dv/dt = (-gL*(v-EL) + gL*DeltaT*exp((v - Vt)/DeltaT) + I - w/S)/Cm : volt (unless refractory) | |
| dw/dt = (a*(v - EL) - w)/tauw : ampere | |
| I = I_inp/S - I_syn: ampere / meter**2 | |
| I_inp = I_app * exp(-r/Lambda) : ampere | |
| r = i*Rad/N : meter | |
| Rad : meter | |
| I_app : ampere | |
| a: siemens | |
| b: ampere | |
| """ | |
| #synapse current equations, same for each neuron group (from Brunel/Wang (ds_Ampa_T/dt) and Gutnisky) | |
| eqs_I_Ampa_T = ''' | |
| I_Ampa_T = g_syn_T * (v - E_Ampa) * s_Ampa_T * 1*ms/(sqrt(Kappa_T)*tau_Ampa): ampere / meter**2 | |
| ds_Ampa_T / dt = - s_Ampa_T / tau_Ampa : 1 | |
| g_syn_T: siemens / meter**2 | |
| Kappa_T :1 | |
| ''' | |
| eqs_I_Ampa_E = ''' | |
| I_Ampa_E = g_syn_E * (v - E_Ampa) * s_Ampa_E * 1*ms/(sqrt(Kappa_E)*tau_Ampa) : ampere / meter**2 | |
| ds_Ampa_E / dt = - s_Ampa_E / tau_Ampa : 1 | |
| g_syn_E: siemens / meter**2 | |
| Kappa_E :1 | |
| ''' | |
| eqs_I_Gaba_I = ''' | |
| I_Gaba_I = g_syn_I * (v - E_Gaba) * s_Gaba_I * 1*ms/(sqrt(Kappa_I)*tau_Gaba): ampere / meter**2 | |
| ds_Gaba_I / dt = - s_Gaba_I / tau_Gaba : 1 | |
| g_syn_I: siemens / meter**2 | |
| Kappa_I :1 | |
| ''' | |
| eqs_I_Gaba_R = ''' | |
| I_Gaba_R = g_syn_R * (v - E_Gaba) * s_Gaba_R * 1*ms/(sqrt(Kappa_R)*tau_Gaba): ampere / meter**2 | |
| ds_Gaba_R / dt = - s_Gaba_R / tau_Gaba : 1 | |
| g_syn_R: siemens / meter**2 | |
| Kappa_R :1 | |
| ''' | |
| #Synaptic Equations, separate for each synapse group | |
| ## Thalamus, Reticular Nrns, L4e, L4i | |
| Thal = NeuronGroup(N_t, | |
| model = eqs_aeIF + | |
| eqs_I_Ampa_E + | |
| eqs_I_Gaba_R + | |
| 'I_syn = (I_Ampa_E + I_Gaba_R): ampere / meter**2', | |
| method = 'euler', threshold = "v>=Vt",reset = 'v=Vr; w+=b', refractory = 2.5*ms) | |
| Thal.a = 0.04 *uS #Destexhe 2009 | |
| Thal.b = 0.0 * nA | |
| Thal.v = Vr | |
| Thal.Rad = R_Th | |
| L4e = NeuronGroup(N_e, | |
| model = eqs_aeIF + | |
| eqs_I_Ampa_T + | |
| eqs_I_Ampa_E + | |
| eqs_I_Gaba_I + | |
| 'I_syn = (I_Ampa_T + I_Ampa_E + I_Gaba_I): ampere / meter**2', | |
| method = 'euler', threshold = "v>=Vt",reset = 'v=Vr; w+=b',refractory = 2.5*ms) | |
| L4e.a = 0.001 * uS | |
| L4e.b = 0.04 * nA | |
| L4e.v = Vr | |
| L4e.Rad = R_Cx | |
| L4i = NeuronGroup(N_i, | |
| model = eqs_aeIF + | |
| eqs_I_Ampa_T + | |
| eqs_I_Ampa_E + | |
| eqs_I_Gaba_I + | |
| 'I_syn = (I_Ampa_T + I_Ampa_E + I_Gaba_I): ampere / meter**2', | |
| method = 'euler',threshold = "v>=Vt",reset = 'v=Vr; w+=b',refractory = 2.5*ms) | |
| L4i.a = 0.001 * uS | |
| L4i.b = 0.0 * nA | |
| L4i.v = Vr | |
| L4i.Rad = R_Cx | |
| Ret = NeuronGroup(N_r, | |
| model = eqs_aeIF + | |
| eqs_I_Ampa_T + | |
| eqs_I_Ampa_E + | |
| eqs_I_Gaba_R + | |
| 'I_syn = (I_Ampa_T + I_Ampa_E + I_Gaba_R): ampere / meter**2', | |
| method = 'euler',threshold = "v>=Vt",reset = 'v=Vr; w+=b',refractory = 2.5*ms) | |
| Ret.a = 0.08 * uS | |
| Ret.b = 0.03 * nA | |
| Ret.v = Vr | |
| Ret.Rad = R_Th | |
| Thal.Kappa_R = Kappa_t_r | |
| Thal.Kappa_E = Kappa_t_e | |
| Thal.g_syn_R = g_T_R | |
| Thal.g_syn_E = g_T_E | |
| Ret.Kappa_T = Kappa_r_t | |
| Ret.Kappa_R = Kappa_r_r | |
| Ret.Kappa_E = Kappa_r_e | |
| Ret.g_syn_T = g_R_T | |
| Ret.g_syn_R = g_R_R | |
| Ret.g_syn_E = g_R_E | |
| L4e.Kappa_T = Kappa_e_t | |
| L4i.Kappa_T = Kappa_i_t | |
| L4e.Kappa_E = Kappa_e_e | |
| L4i.Kappa_E = Kappa_i_e | |
| L4e.Kappa_I = Kappa_e_i | |
| L4i.Kappa_I = Kappa_i_i | |
| L4e.g_syn_T = g_E_T | |
| L4i.g_syn_T = g_I_T | |
| L4e.g_syn_E = g_E_E | |
| L4i.g_syn_E = g_I_E | |
| L4e.g_syn_I = g_E_I | |
| L4i.g_syn_I = g_I_I | |
| P_t_exc = PoissonInput(Thal, 'v', Po_N_exc, Po_Rate, weight=Po_Amp_T) | |
| P_i_exc = PoissonInput(L4i, 'v', Po_N_exc, Po_Rate, weight=Po_Amp_I) | |
| P_e_exc = PoissonInput(L4e, 'v', Po_N_exc, Po_Rate, weight=Po_Amp_E) | |
| P_r_exc = PoissonInput(Ret, 'v', Po_N_exc, Po_Rate, weight=Po_Amp_R) | |
| P_t_inh = PoissonInput(Thal, 'v', Po_N_inh, Po_Rate, weight=-Po_Amp_T) | |
| P_i_inh = PoissonInput(L4i, 'v', Po_N_inh, Po_Rate, weight=-Po_Amp_I) | |
| P_e_inh = PoissonInput(L4e, 'v', Po_N_inh, Po_Rate, weight=-Po_Amp_E) | |
| P_r_inh = PoissonInput(Ret, 'v', Po_N_inh, Po_Rate, weight=-Po_Amp_R) | |
| #%% Connectivity and synapses | |
| #read as: connect thalamus to l4e (no connectivity spec = fully connected) | |
| S_e_t = Synapses(Thal, L4e, on_pre='s_Ampa_T+=1', delay = 1*ms) | |
| S_i_t = Synapses(Thal, L4i, on_pre='s_Ampa_T+=1', delay = 1*ms) | |
| S_e_e = Synapses(L4e, L4e, on_pre='s_Ampa_E+=1', delay = 1*ms) | |
| S_i_e = Synapses(L4e, L4i, on_pre='s_Ampa_E+=1', delay = 1*ms) | |
| S_e_i = Synapses(L4i, L4e, on_pre='s_Gaba_I+=1', delay = .85*ms) | |
| S_i_i = Synapses(L4i, L4i, on_pre='s_Gaba_I+=1', delay = 0.5*ms) | |
| S_t_e = Synapses(L4e, Thal, on_pre='s_Ampa_E+=1', delay = 1*ms) | |
| S_r_e = Synapses(L4e, Ret, on_pre='s_Ampa_E+=1', delay = 1*ms) | |
| S_r_t = Synapses(Thal, Ret, on_pre='s_Ampa_T+=1', delay = 1*ms) | |
| S_t_r = Synapses(Ret, Thal, on_pre='s_Gaba_R+=1', delay = 0.5*ms) | |
| S_r_r = Synapses(Ret, Ret, on_pre='s_Gaba_R+=1', delay = 0.5*ms) | |
| S_e_t.connect(p=Kappa_e_t/N_t) # 200 T -> 1600 E (25%) | |
| S_i_t.connect(p=Kappa_i_t/N_t) | |
| S_e_e.connect(p=Kappa_e_e/N_e) | |
| S_i_e.connect(p=Kappa_i_e/N_e) | |
| S_e_i.connect(p=Kappa_e_i/N_i) | |
| S_i_i.connect(p=Kappa_i_i/N_i) | |
| S_t_e.connect(p=Kappa_t_e/N_e) # 1600 E -> 200 T (2%) | |
| S_r_e.connect(p=Kappa_r_e/N_e) | |
| S_r_t.connect(p=Kappa_r_t/N_t) | |
| S_t_r.connect(p=Kappa_t_r/N_r) | |
| S_r_r.connect(p=Kappa_r_r/N_r) | |
| # network definition | |
| Thal.I_app = 0*pA | |
| L4e.I_app = 0*pA | |
| L4i.I_app = 0*pA | |
| Ret.I_app = 0*pA | |
| Spike_t = SpikeMonitor(Thal) | |
| Spike_e = SpikeMonitor(L4e) | |
| Spike_i = SpikeMonitor(L4i) | |
| Spike_r = SpikeMonitor(Ret) | |
| ratE_t = PopulationRateMonitor(Thal) | |
| ratE_e = PopulationRateMonitor(L4e) | |
| ratE_i = PopulationRateMonitor(L4i) | |
| ratE_r = PopulationRateMonitor(Ret) | |
| run(Init_time, report='stdout') | |
| fig, axs = plt.subplots(1, 2, figsize=(10, 4), sharex="col", sharey="row") | |
| for i, (rate_mon, title) in enumerate(zip([ratE_t, ratE_r, ratE_e], ['Thalamus', 'Ret', 'L4e'])): | |
| axs[0].plot(rate_mon.t/second, rate_mon.rate/Hz, label=title) | |
| axs[0].set_xlim(48.8, 50) | |
| axs[0].legend() | |
| store() | |
| # a few steps dedicated to establishing directories for output data | |
| i, s = 0, 0 | |
| restore() | |
| del Spike_t, Spike_e, Spike_i, Spike_r | |
| del ratE_t, ratE_e, ratE_i, ratE_r | |
| Spike_t = SpikeMonitor(Thal) | |
| Spike_e = SpikeMonitor(L4e) | |
| Spike_i = SpikeMonitor(L4i) | |
| Spike_r = SpikeMonitor(Ret) | |
| ratE_t = PopulationRateMonitor(Thal) | |
| ratE_e = PopulationRateMonitor(L4e) | |
| ratE_i = PopulationRateMonitor(L4i) | |
| ratE_r = PopulationRateMonitor(Ret) | |
| setParam('Istim', str(i) + '*pA', params) | |
| setParam('StimCx', str(s), params) | |
| Thal.I_app = Istim * (1-StimCx) | |
| L4e.I_app = Istim * StimCx | |
| L4i.I_app = Istim * StimCx | |
| run(run_time, report='stdout') | |
| for i, (rate_mon, title) in enumerate(zip([ratE_t, ratE_r, ratE_e], ['Thalamus', 'Ret', 'L4e'])): | |
| axs[1].plot(rate_mon.t/second, rate_mon.rate/Hz, label=title) | |
| axs[1].set_xlim(50, 51.2) | |
| axs[1].legend() | |
| plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment