.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/gif_pop_psc_exp.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_gif_pop_psc_exp.py: Population rate model of generalized integrate-and-fire neurons --------------------------------------------------------------- .. only:: html ---- Run this example as a Jupyter notebook: .. card:: :width: 25% :margin: 2 :text-align: center :link: https://lab.ebrains.eu/hub/user-redirect/git-pull?repo=https%3A%2F%2Fgithub.com%2Fnest%2Fnest-simulator-examples&urlpath=lab%2Ftree%2Fnest-simulator-examples%2Fnotebooks%2Fnotebooks%2Fgif_pop_psc_exp.ipynb&branch=main :link-alt: JupyterHub service .. image:: https://nest-simulator.org/TryItOnEBRAINS.png .. grid:: 1 1 1 1 :padding: 0 0 2 0 .. grid-item:: :class: sd-text-muted :margin: 0 0 3 0 :padding: 0 0 3 0 :columns: 4 See :ref:`our guide ` for more information and troubleshooting. ---- This script simulates a finite network of generalized integrate-and-fire (GIF) neurons directly on the mesoscopic population level using the effective stochastic population rate dynamics derived in the paper [1]_. The stochastic population dynamics is implemented in the NEST model gif_pop_psc_exp. We demonstrate this model using the example of a Brunel network of two coupled populations, one excitatory and one inhibitory population. Note that the population model represents the mesoscopic level description of the corresponding microscopic network based on the NEST model ``gif_psc_exp``. References ~~~~~~~~~~ .. [1] Schwalger T, Degert M, Gerstner W (2017). Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. PLoS Comput Biol. https://doi.org/10.1371/journal.pcbi.1005507 .. GENERATED FROM PYTHON SOURCE LINES 48-49 Import necessary modules. .. GENERATED FROM PYTHON SOURCE LINES 49-54 .. code-block:: Python import matplotlib.pyplot as plt import nest import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 55-56 We first set the parameters of the microscopic model: .. GENERATED FROM PYTHON SOURCE LINES 56-111 .. code-block:: Python # All times given in milliseconds dt = 0.5 dt_rec = 1.0 # Simulation time t_end = 2000.0 # Parameters size = 200 N = np.array([4, 1]) * size M = len(N) # number of populations # neuronal parameters t_ref = 4.0 * np.ones(M) # absolute refractory period tau_m = 20 * np.ones(M) # membrane time constant mu = 24.0 * np.ones(M) # constant base current mu=R*(I0+Vrest) c = 10.0 * np.ones(M) # base rate of exponential link function Delta_u = 2.5 * np.ones(M) # softness of exponential link function V_reset = 0.0 * np.ones(M) # Reset potential V_th = 15.0 * np.ones(M) # baseline threshold (non-accumulating part) tau_sfa_exc = [100.0, 1000.0] # adaptation time constants of excitatory neurons tau_sfa_inh = [100.0, 1000.0] # adaptation time constants of inhibitory neurons J_sfa_exc = [1000.0, 1000.0] # size of feedback kernel theta # (= area under exponential) in mV*ms J_sfa_inh = [1000.0, 1000.0] # in mV*ms tau_theta = np.array([tau_sfa_exc, tau_sfa_inh]) J_theta = np.array([J_sfa_exc, J_sfa_inh]) # connectivity J = 0.3 # excitatory synaptic weight in mV if number of input connections # is C0 (see below) g = 5.0 # inhibition-to-excitation ratio pconn = 0.2 * np.ones((M, M)) delay = 1.0 * np.ones((M, M)) C0 = np.array([[800, 200], [800, 200]]) * 0.2 # constant reference matrix C = np.vstack((N, N)) * pconn # numbers of input connections # final synaptic weights scaling as 1/C J_syn = np.array([[J, -g * J], [J, -g * J]]) * C0 / C taus1_ = [3.0, 6.0] # time constants of exc./inh. postsynaptic currents (PSC's) taus1 = np.array([taus1_ for k in range(M)]) # step current input step = [[20.0], [20.0]] # jump size of mu in mV tstep = np.array([[1500.0], [1500.0]]) # times of jumps # synaptic time constants of excitatory and inhibitory connections tau_ex = 3.0 # in ms tau_in = 6.0 # in ms .. GENERATED FROM PYTHON SOURCE LINES 112-119 Simulation on the mesoscopic level ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ To directly simulate the mesoscopic population activities (i.e. generating the activity of a finite-size population without simulating single neurons), we can build the populations using the NEST model ``gif_pop_psc_exp``: .. GENERATED FROM PYTHON SOURCE LINES 119-167 .. code-block:: Python nest.set_verbosity("M_WARNING") nest.ResetKernel() nest.resolution = dt nest.print_time = True nest.local_num_threads = 1 t0 = nest.biological_time nest_pops = nest.Create("gif_pop_psc_exp", M) C_m = 250.0 # irrelevant value for membrane capacity, cancels out in simulation g_L = C_m / tau_m params = [ { "C_m": C_m, "I_e": mu[i] * g_L[i], "lambda_0": c[i], # in Hz! "Delta_V": Delta_u[i], "tau_m": tau_m[i], "tau_sfa": tau_theta[i], "q_sfa": J_theta[i] / tau_theta[i], # [J_theta]= mV*ms -> [q_sfa]=mV "V_T_star": V_th[i], "V_reset": V_reset[i], "len_kernel": -1, # -1 triggers automatic history size "N": N[i], "t_ref": t_ref[i], "tau_syn_ex": max([tau_ex, dt]), "tau_syn_in": max([tau_in, dt]), "E_L": 0.0, } for i in range(M) ] nest_pops.set(params) # connect the populations g_syn = np.ones_like(J_syn) # synaptic conductance g_syn[:, 0] = C_m / tau_ex g_syn[:, 1] = C_m / tau_in for i in range(M): for j in range(M): nest.Connect( nest_pops[j], nest_pops[i], syn_spec={"weight": J_syn[i, j] * g_syn[i, j] * pconn[i, j], "delay": delay[i, j]}, ) .. GENERATED FROM PYTHON SOURCE LINES 168-170 To record the instantaneous population rate `Abar(t)` we use a multimeter, and to get the population activity `A_N(t)` we use spike recorder: .. GENERATED FROM PYTHON SOURCE LINES 170-183 .. code-block:: Python # monitor the output using a multimeter, this only records with dt_rec! nest_mm = nest.Create("multimeter") nest_mm.set(record_from=["n_events", "mean"], interval=dt_rec) nest.Connect(nest_mm, nest_pops) # monitor the output using a spike recorder nest_sr = [] for i in range(M): nest_sr.append(nest.Create("spike_recorder")) nest_sr[i].time_in_steps = True nest.Connect(nest_pops[i], nest_sr[i], syn_spec={"weight": 1.0, "delay": dt}) .. GENERATED FROM PYTHON SOURCE LINES 184-186 All neurons in a given population will be stimulated with a step input current: .. GENERATED FROM PYTHON SOURCE LINES 186-199 .. code-block:: Python # set initial value (at t0+dt) of step current generator to zero tstep = np.hstack((dt * np.ones((M, 1)), tstep)) step = np.hstack((np.zeros((M, 1)), step)) # create the step current devices nest_stepcurrent = nest.Create("step_current_generator", M) # set the parameters for the step currents for i in range(M): nest_stepcurrent[i].set(amplitude_times=tstep[i] + t0, amplitude_values=step[i] * g_L[i], origin=t0, stop=t_end) pop_ = nest_pops[i] nest.Connect(nest_stepcurrent[i], pop_, syn_spec={"weight": 1.0, "delay": dt}) .. GENERATED FROM PYTHON SOURCE LINES 200-201 We can now start the simulation: .. GENERATED FROM PYTHON SOURCE LINES 201-223 .. code-block:: Python nest.rng_seed = 1 t = np.arange(0.0, t_end, dt_rec) A_N = np.ones((t.size, M)) * np.nan Abar = np.ones_like(A_N) * np.nan # simulate 1 step longer to make sure all t are simulated nest.Simulate(t_end + dt) data_mm = nest_mm.events for i, nest_i in enumerate(nest_pops): a_i = data_mm["mean"][data_mm["senders"] == nest_i.global_id] a = a_i / N[i] / dt min_len = np.min([len(a), len(Abar)]) Abar[:min_len, i] = a[:min_len] data_sr = nest_sr[i].get("events", "times") data_sr = data_sr * dt - t0 bins = np.concatenate((t, np.array([t[-1] + dt_rec]))) A = np.histogram(data_sr, bins=bins)[0] / float(N[i]) / dt_rec A_N[:, i] = A .. GENERATED FROM PYTHON SOURCE LINES 224-225 and plot the activity: .. GENERATED FROM PYTHON SOURCE LINES 225-237 .. code-block:: Python plt.figure(1) plt.clf() plt.subplot(2, 1, 1) plt.plot(t, A_N * 1000) # plot population activities (in Hz) plt.ylabel(r"$A_N$ [Hz]") plt.title("Population activities (mesoscopic sim.)") plt.subplot(2, 1, 2) plt.plot(t, Abar * 1000) # plot instantaneous population rates (in Hz) plt.ylabel(r"$\bar A$ [Hz]") plt.xlabel("time [ms]") .. GENERATED FROM PYTHON SOURCE LINES 238-246 Microscopic ("direct") simulation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ As mentioned above, the population model ``gif_pop_psc_exp`` directly simulates the mesoscopic population activities, i.e. without the need to simulate single neurons. On the other hand, if we want to know single neuron activities, we must simulate on the microscopic level. This is possible by building a corresponding network of ``gif_psc_exp`` neuron models: .. GENERATED FROM PYTHON SOURCE LINES 246-287 .. code-block:: Python nest.ResetKernel() nest.resolution = dt nest.print_time = True nest.local_num_threads = 1 t0 = nest.biological_time nest_pops = [] for k in range(M): nest_pops.append(nest.Create("gif_psc_exp", N[k])) # set single neuron properties for i in range(M): nest_pops[i].set( C_m=C_m, I_e=mu[i] * g_L[i], lambda_0=c[i], Delta_V=Delta_u[i], g_L=g_L[i], tau_sfa=tau_theta[i], q_sfa=J_theta[i] / tau_theta[i], V_T_star=V_th[i], V_reset=V_reset[i], t_ref=t_ref[i], tau_syn_ex=max([tau_ex, dt]), tau_syn_in=max([tau_in, dt]), E_L=0.0, V_m=0.0, ) # connect the populations for i, nest_i in enumerate(nest_pops): for j, nest_j in enumerate(nest_pops): if np.allclose(pconn[i, j], 1.0): conn_spec = {"rule": "all_to_all"} else: conn_spec = {"rule": "fixed_indegree", "indegree": int(pconn[i, j] * N[j])} nest.Connect(nest_j, nest_i, conn_spec, syn_spec={"weight": J_syn[i, j] * g_syn[i, j], "delay": delay[i, j]}) .. GENERATED FROM PYTHON SOURCE LINES 288-291 We want to record all spikes of each population in order to compute the mesoscopic population activities `A_N(t)` from the microscopic simulation. We also record the membrane potentials of five example neurons: .. GENERATED FROM PYTHON SOURCE LINES 291-309 .. code-block:: Python # monitor the output using a multimeter and a spike recorder nest_sr = [] for i, nest_i in enumerate(nest_pops): nest_sr.append(nest.Create("spike_recorder")) nest_sr[i].time_in_steps = True # record all spikes from population to compute population activity nest.Connect(nest_i, nest_sr[i], syn_spec={"weight": 1.0, "delay": dt}) Nrecord = [5, 0] # for each population "i" the first Nrecord[i] neurons are recorded nest_mm_Vm = [] for i, nest_i in enumerate(nest_pops): nest_mm_Vm.append(nest.Create("multimeter")) nest_mm_Vm[i].set(record_from=["V_m"], interval=dt_rec) if Nrecord[i] != 0: nest.Connect(nest_mm_Vm[i], nest_i[: Nrecord[i]], syn_spec={"weight": 1.0, "delay": dt}) .. GENERATED FROM PYTHON SOURCE LINES 310-313 As before, all neurons in a given population will be stimulated with a step input current. The following code block is identical to the one for the mesoscopic simulation above: .. GENERATED FROM PYTHON SOURCE LINES 313-324 .. code-block:: Python # create the step current devices if they do not exist already nest_stepcurrent = nest.Create("step_current_generator", M) # set the parameters for the step currents for i in range(M): nest_stepcurrent[i].set(amplitude_times=tstep[i] + t0, amplitude_values=step[i] * g_L[i], origin=t0, stop=t_end) nest_stepcurrent[i].set(amplitude_times=tstep[i] + t0, amplitude_values=step[i] * g_L[i], origin=t0, stop=t_end) # optionally a stopping time may be added by: 'stop': sim_T + t0 pop_ = nest_pops[i] nest.Connect(nest_stepcurrent[i], pop_, syn_spec={"weight": 1.0, "delay": dt}) .. GENERATED FROM PYTHON SOURCE LINES 325-326 We can now start the microscopic simulation: .. GENERATED FROM PYTHON SOURCE LINES 326-335 .. code-block:: Python nest.rng_seed = 1 t = np.arange(0.0, t_end, dt_rec) A_N = np.ones((t.size, M)) * np.nan # simulate 1 step longer to make sure all t are simulated nest.Simulate(t_end + dt) .. GENERATED FROM PYTHON SOURCE LINES 336-338 Let's retrieve the data of the spike recorder and plot the activity of the excitatory population (in Hz): .. GENERATED FROM PYTHON SOURCE LINES 338-352 .. code-block:: Python for i in range(len(nest_pops)): data_sr = nest_sr[i].get("events", "times") * dt - t0 bins = np.concatenate((t, np.array([t[-1] + dt_rec]))) A = np.histogram(data_sr, bins=bins)[0] / float(N[i]) / dt_rec A_N[:, i] = A * 1000 # in Hz t = np.arange(dt, t_end + dt, dt_rec) plt.figure(2) plt.plot(t, A_N[:, 0]) plt.xlabel("time [ms]") plt.ylabel("population activity [Hz]") plt.title("Population activities (microscopic sim.)") .. GENERATED FROM PYTHON SOURCE LINES 353-358 This should look similar to the population activity obtained from the mesoscopic simulation based on the NEST model ``gif_pop_psc_exp`` (cf. figure 1). Now we retrieve the data of the multimeter, which allows us to look at the membrane potentials of single neurons. Here we plot the voltage traces (in mV) of five example neurons: .. GENERATED FROM PYTHON SOURCE LINES 358-376 .. code-block:: Python voltage = [] for i in range(M): if Nrecord[i] > 0: senders = nest_mm_Vm[i].get("events", "senders") v = nest_mm_Vm[i].get("events", "V_m") voltage.append(np.array([v[np.where(senders == j)] for j in set(senders)])) else: voltage.append(np.array([])) f, axarr = plt.subplots(Nrecord[0], sharex=True) for i in range(Nrecord[0]): axarr[i].plot(voltage[0][i]) axarr[i].set_yticks((0, 15, 30)) axarr[i].set_xlabel("time [ms]") axarr[2].set_ylabel("membrane potential [mV]") axarr[0].set_title("5 example GIF neurons (microscopic sim.)") .. GENERATED FROM PYTHON SOURCE LINES 377-379 Note that this plots only the subthreshold membrane potentials but not the spikes (as with every leaky integrate-and-fire model). .. GENERATED FROM PYTHON SOURCE LINES 379-381 .. code-block:: Python plt.show() .. _sphx_glr_download_auto_examples_gif_pop_psc_exp.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: gif_pop_psc_exp.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: gif_pop_psc_exp.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_