.. currentmodule:: brian2 .. Brunel_Wang_2001: Example: Brunel_Wang_2001 ========================= .. only:: html .. |launchbinder| image:: http://mybinder.org/badge.svg .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/frompapers/Brunel_Wang_2001.ipynb .. note:: You can launch an interactive, editable version of this example without installing any local files using the Binder service (although note that at some times this may be slow or fail to open): |launchbinder|_ Sample-specific persistent activity ----------------------------------- Five subpopulations with three selective and one reset stimuli example. Analog to figure 6b in the paper. BRUNEL, Nicolas et WANG, Xiao-Jing. Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition. Journal of computational neuroscience, 2001, vol. 11, no 1, p. 63-85. :: from brian2 import * # populations N = 1000 N_E = int(N * 0.8) # pyramidal neurons N_I = int(N * 0.2) # interneurons # voltage V_L = -70. * mV V_thr = -50. * mV V_reset = -55. * mV V_E = 0. * mV V_I = -70. * mV # membrane capacitance C_m_E = 0.5 * nF C_m_I = 0.2 * nF # membrane leak g_m_E = 25. * nS g_m_I = 20. * nS # refractory period tau_rp_E = 2. * ms tau_rp_I = 1. * ms # external stimuli rate = 3 * Hz C_ext = 800 # synapses C_E = N_E C_I = N_I # AMPA (excitatory) g_AMPA_ext_E = 2.08 * nS g_AMPA_rec_E = 0.104 * nS * 800. / N_E g_AMPA_ext_I = 1.62 * nS g_AMPA_rec_I = 0.081 * nS * 800. / N_E tau_AMPA = 2. * ms # NMDA (excitatory) g_NMDA_E = 0.327 * nS * 800. / N_E g_NMDA_I = 0.258 * nS * 800. / N_E tau_NMDA_rise = 2. * ms tau_NMDA_decay = 100. * ms alpha = 0.5 / ms Mg2 = 1. # GABAergic (inhibitory) g_GABA_E = 1.25 * nS * 200. / N_I g_GABA_I = 0.973 * nS * 200. / N_I tau_GABA = 10. * ms # subpopulations f = 0.1 p = 5 N_sub = int(N_E * f) N_non = int(N_E * (1. - f * p)) w_plus = 2.1 w_minus = 1. - f * (w_plus - 1.) / (1. - f) # modeling eqs_E = ''' dv / dt = (- g_m_E * (v - V_L) - I_syn) / C_m_E : volt (unless refractory) I_syn = I_AMPA_ext + I_AMPA_rec + I_NMDA_rec + I_GABA_rec : amp I_AMPA_ext = g_AMPA_ext_E * (v - V_E) * s_AMPA_ext : amp I_AMPA_rec = g_AMPA_rec_E * (v - V_E) * 1 * s_AMPA : amp ds_AMPA_ext / dt = - s_AMPA_ext / tau_AMPA : 1 ds_AMPA / dt = - s_AMPA / tau_AMPA : 1 I_NMDA_rec = g_NMDA_E * (v - V_E) / (1 + Mg2 * exp(-0.062 * v / mV) / 3.57) * s_NMDA_tot : amp s_NMDA_tot : 1 I_GABA_rec = g_GABA_E * (v - V_I) * s_GABA : amp ds_GABA / dt = - s_GABA / tau_GABA : 1 ''' eqs_I = ''' dv / dt = (- g_m_I * (v - V_L) - I_syn) / C_m_I : volt (unless refractory) I_syn = I_AMPA_ext + I_AMPA_rec + I_NMDA_rec + I_GABA_rec : amp I_AMPA_ext = g_AMPA_ext_I * (v - V_E) * s_AMPA_ext : amp I_AMPA_rec = g_AMPA_rec_I * (v - V_E) * 1 * s_AMPA : amp ds_AMPA_ext / dt = - s_AMPA_ext / tau_AMPA : 1 ds_AMPA / dt = - s_AMPA / tau_AMPA : 1 I_NMDA_rec = g_NMDA_I * (v - V_E) / (1 + Mg2 * exp(-0.062 * v / mV) / 3.57) * s_NMDA_tot : amp s_NMDA_tot : 1 I_GABA_rec = g_GABA_I * (v - V_I) * s_GABA : amp ds_GABA / dt = - s_GABA / tau_GABA : 1 ''' P_E = NeuronGroup(N_E, eqs_E, threshold='v > V_thr', reset='v = V_reset', refractory=tau_rp_E, method='euler') P_E.v = V_L P_I = NeuronGroup(N_I, eqs_I, threshold='v > V_thr', reset='v = V_reset', refractory=tau_rp_I, method='euler') P_I.v = V_L eqs_glut = ''' s_NMDA_tot_post = w * s_NMDA : 1 (summed) ds_NMDA / dt = - s_NMDA / tau_NMDA_decay + alpha * x * (1 - s_NMDA) : 1 (clock-driven) dx / dt = - x / tau_NMDA_rise : 1 (clock-driven) w : 1 ''' eqs_pre_glut = ''' s_AMPA += w x += 1 ''' eqs_pre_gaba = ''' s_GABA += 1 ''' # E to E C_E_E = Synapses(P_E, P_E, model=eqs_glut, on_pre=eqs_pre_glut, method='euler') C_E_E.connect('i != j') C_E_E.w[:] = 1 for pi in range(N_non, N_non + p * N_sub, N_sub): # internal other subpopulation to current nonselective C_E_E.w[C_E_E.indices[:, pi:pi + N_sub]] = w_minus # internal current subpopulation to current subpopulation C_E_E.w[C_E_E.indices[pi:pi + N_sub, pi:pi + N_sub]] = w_plus # E to I C_E_I = Synapses(P_E, P_I, model=eqs_glut, on_pre=eqs_pre_glut, method='euler') C_E_I.connect() C_E_I.w[:] = 1 # I to I C_I_I = Synapses(P_I, P_I, on_pre=eqs_pre_gaba, method='euler') C_I_I.connect('i != j') # I to E C_I_E = Synapses(P_I, P_E, on_pre=eqs_pre_gaba, method='euler') C_I_E.connect() # external noise C_P_E = PoissonInput(P_E, 's_AMPA_ext', C_ext, rate, '1') C_P_I = PoissonInput(P_I, 's_AMPA_ext', C_ext, rate, '1') # at 1s, select population 1 C_selection = int(f * C_ext) rate_selection = 25 * Hz stimuli1 = TimedArray(np.r_[np.zeros(40), np.ones(2), np.zeros(100)], dt=25 * ms) input1 = PoissonInput(P_E[N_non:N_non + N_sub], 's_AMPA_ext', C_selection, rate_selection, 'stimuli1(t)') # at 2s, select population 2 stimuli2 = TimedArray(np.r_[np.zeros(80), np.ones(2), np.zeros(100)], dt=25 * ms) input2 = PoissonInput(P_E[N_non + N_sub:N_non + 2 * N_sub], 's_AMPA_ext', C_selection, rate_selection, 'stimuli2(t)') # at 4s, reset selection stimuli_reset = TimedArray(np.r_[np.zeros(120), np.ones(2), np.zeros(100)], dt=25 * ms) input_reset_I = PoissonInput(P_E, 's_AMPA_ext', C_ext, rate_selection, 'stimuli_reset(t)') input_reset_E = PoissonInput(P_I, 's_AMPA_ext', C_ext, rate_selection, 'stimuli_reset(t)') # monitors N_activity_plot = 15 sp_E_sels = [SpikeMonitor(P_E[pi:pi + N_activity_plot]) for pi in range(N_non, N_non + p * N_sub, N_sub)] sp_E = SpikeMonitor(P_E[:N_activity_plot]) sp_I = SpikeMonitor(P_I[:N_activity_plot]) r_E_sels = [PopulationRateMonitor(P_E[pi:pi + N_sub]) for pi in range(N_non, N_non + p * N_sub, N_sub)] r_E = PopulationRateMonitor(P_E[:N_non]) r_I = PopulationRateMonitor(P_I) # simulate, can be long >120s net = Network(collect()) net.add(sp_E_sels) net.add(r_E_sels) net.run(4 * second, report='stdout') # plotting title('Population rates') xlabel('ms') ylabel('Hz') plot(r_E.t / ms, r_E.smooth_rate(width=25 * ms) / Hz, label='nonselective') plot(r_I.t / ms, r_I.smooth_rate(width=25 * ms) / Hz, label='inhibitory') for i, r_E_sel in enumerate(r_E_sels[::-1]): plot(r_E_sel.t / ms, r_E_sel.smooth_rate(width=25 * ms) / Hz, label=f"selective {p - i}") legend() figure() title(f"Population activities ({N_activity_plot} neurons/pop)") xlabel('ms') yticks([]) plot(sp_E.t / ms, sp_E.i + (p + 1) * N_activity_plot, '.', markersize=2, label="nonselective") plot(sp_I.t / ms, sp_I.i + p * N_activity_plot, '.', markersize=2, label="inhibitory") for i, sp_E_sel in enumerate(sp_E_sels[::-1]): plot(sp_E_sel.t / ms, sp_E_sel.i + (p - i - 1) * N_activity_plot, '.', markersize=2, label=f"selective {p - i}") legend() show() .. image:: ../resources/examples_images/frompapers.Brunel_Wang_2001.1.png .. image:: ../resources/examples_images/frompapers.Brunel_Wang_2001.2.png