10. Network of LIF neurons (Brunel)

In this exercise we study a well known network of sparsely connected Leaky-Integrate-And-Fire neurons (Brunel, 2000).

Book chapters

The Brunel model is introduced in Chapter 13 Section 4.2. The network structure is shown in Figure 13.6b. Read the section Synchrony, oscillations, and irregularity and have a look at Figure 13.7. For this exercise, you can skip the explanations related to the Fokker-Planck equation.

Python classes

The module brunel_model.LIF_spiking_network implements a parametrized network. The figure below shows the simulation result using the default configuration.

../_images/Brunel_Spiking_LIF.png

Simulation result. Top: raster plot of 150 randomly selected neurons. Three spike trains are visually highlighted. Middle: time evolution of the population activity \(A(t)\). Bottom: Membrane voltage of three neurons. The red color in the top and bottom panels identifies the same neuron.

To get started, call the function brunel_model.LIF_spiking_network.getting_started() or copy the following code into a Jupyter notebook.

%matplotlib inline
from neurodynex3.brunel_model import LIF_spiking_network
from neurodynex3.tools import plot_tools
import brian2 as b2

rate_monitor, spike_monitor, voltage_monitor, monitored_spike_idx = LIF_spiking_network.simulate_brunel_network(sim_time=250. * b2.ms)
plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min=0.*b2.ms)

Note that you can change all parameters of the neuron by using the named parameters of the function simulate_brunel_network(). If you do not specify any parameter, the default values are used (see next code block). You can access these variables in your code by prefixing them with the module name (for example LIF_spiking_network.POISSON_INPUT_RATE).

# Default parameters of a single LIF neuron:
V_REST = 0. * b2.mV
V_RESET = +10. * b2.mV
FIRING_THRESHOLD = +20. * b2.mV
MEMBRANE_TIME_SCALE = 20. * b2.ms
ABSOLUTE_REFRACTORY_PERIOD = 2.0 * b2.ms

# Default parameters of the network
SYNAPTIC_WEIGHT_W0 = 0.1 * b2.mV  # note: w_ee=w_ie = w0 and = w_ei=w_ii = -g*w0
RELATIVE_INHIBITORY_STRENGTH_G = 4.  # balanced
CONNECTION_PROBABILITY_EPSILON = 0.1
SYNAPTIC_DELAY = 1.5 * b2.ms
POISSON_INPUT_RATE = 12. * b2.Hz
N_POISSON_INPUT = 1000

10.1. Exercise: model parameters and threshold rate

In the first exercise, we get familiar with the model and parameters. Make sure you have read the book chapter . Then have a look at the documentation of simulate_brunel_network(). Note that in our implementation, the number of excitatory presynaptic poisson neurons (input from the external population) is a parameter \(N_extern\) and thus independent of \(C_E\).

10.1.1. Question:

  • Run the simulation with the default parameters (see code block above). In that default configuration, what are the values of the variables \(N_E\), \(N_I\), \(C_E\), \(C_I\), \(w_{EE}\), \(w_{EI}\), \(w_{IE}\), and \(w_{II}\)? The variables are described in the book and in Fig. 13.6.
  • What are the units of the weights \(w\)?
  • The frequency \(\nu_{threshold}\) is is the poisson rate of the external population sufficient to drive the neurons in the network to the firing threshold. Using Eq. (1), compute \(\nu_{threshold}\). You can do this in Python, e.g. use LIF_spiking_network.FIRING_THRESHOLD for \(u_{thr}\), etc.
  • Refering to Fig. 13.7, left panel, what is the meaning of the value 1 on the y-axis (Input). What is the horizontal dashed line designating? How is it related to \(u_{thr}\)?
  • Run a simulation for 500ms. Set poisson_input_rate to \(\nu_{threshold}\). Plot the network activity in the time interval [0ms, 500ms]. Is the network quiet (Q)?
  • During the simulation time, what is the average firing rate of a single neuron? You can access the total number of spikes from the SpikeMonitor object in Brian2 with the command spike_monitor.num_spikes and the number of neurons in the network with spike_monitor.source.N.
(1)\[ \nu_{threshold} = \frac{u_{thr}}{N_{extern} w_{0} \tau_m}\]

10.2. Exercise: Population activity

The network of spiking LIF-neurons shows characteristic population activities. In this exercise we investigate the asynchronous irregular (AI), synchronous regular (SR), fast synchronous irregular (SI-fast) and slow synchronous irregular (SI-slow) activity types.

10.2.1. Question: Network states

  • The function simulate_brunel_network() gives you three options to vary the input strength (y-axis in Figure 13.7a). Which options do you have?
  • Which parameter of the function simulate_brunel_network() lets you change the relative strength of inhibition (the x-axis in Figure 13.7, a)?
  • Define a network of 6000 excitatory and 1500 inhibitory neurons. Find the appropriate parameters and simulate the network in the regimes AI, SR, SI-fast and SI-slow. For each of the four configurations, plot the network activity and compute the average firing rate. Run each simulation for at least 1000ms and plot two figures for each simulation: one showing the complete simulation time and one showing only the last ~50ms.
  • What is the population activity \(A(t)\) in each of the four conditions (in Hz, averaged over the last 200ms of your simulation)?

10.2.2. Question: Interspike interval (ISI) and Coefficient of Variation (CV)

Before answering the questions, make sure you understand the notions ISI and CV. If necessary, read Chapter 7.3.1 .

  • What is the CV of a Poisson neuron?

  • From the four figures plotted in the previous question, qualitatively interpret the spike trains and the population activity in each of the four regimes:

    • What is the mean firing rate of a single neuron (only a rough estimate).
    • Sketch the ISI histogram. (Is it peaked or broad? Where’s the maximum?)
    • Estimate the CV. (Is it \(<1\), \(\ll 1\), \(=1\), \(>1\)?)
  • Validate your estimates using the functions spike_tools.get_spike_train_stats() and plot_tools.plot_ISI_distribution(). Use the code block provided here.

  • Make sure you understand the code block. Why is the function spike_tools.get_spike_train_stats called with the parameter window_t_min=100.*b2.ms?

%matplotlib inline
from neurodynex3.brunel_model import LIF_spiking_network
from neurodynex3.tools import plot_tools, spike_tools
import brian2 as b2

poisson_rate = XXXX *b2.Hz
g = XXXX
CE = XXXX
simtime = XXXX *b2.ms

rate_monitor, spike_monitor, voltage_monitor, monitored_spike_idx = LIF_spiking_network.simulate_brunel_network(N_Excit=CE, poisson_input_rate=poisson_rate, g=g, sim_time=simtime)
plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min = 0*b2.ms)
plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min = simtime - XXXX *b2.ms)
spike_stats = spike_tools.get_spike_train_stats(spike_monitor, window_t_min= 100 *b2.ms)
plot_tools.plot_ISI_distribution(spike_stats, hist_nr_bins=100, xlim_max_ISI= XXXX *b2.ms)
  • In the Synchronous Regular (SR) state, what is the dominant frequency of the population activity \(A(t)\)? Compare this frequency to the firing frequency of a single neuron. You can do this “visually” using the plots created by plot_tools.plot_network_activity() or by solving the bonus exercise below.

10.3. Exercise: Emergence of Synchronization

The different regimes emerge from the recurrence and the relative strength of inhibition \(g\). In the absence of recurrent excitatory connections, the network would approach a constant mean activity \(A(t)\).

10.3.1. Question:

  • Simulate a network of 6000 excitatory and 1500 inhibitory neurons. Set the following parameters: poisson_rate=14*b2.Hz, g=2.5. In which state is this network?
  • What would the population activity be if we would have only external input? We can simulate this. Run a simulation of the same network, but disable the recurrent excitatory connections: simulate_brunel_network(...,w0=0.*b2.mV, w_external=LIF_spiking_network.SYNAPTIC_WEIGHT_W0).
  • Explain why the non-recurrent network shows a strong synchronization in the beginning and why this synchronization fades out.
  • The non-recurrent network is strongly synchronized in the beginning. Is the connected network simply “locked” to this initial synchronization? You can falsify this hypothesis by initializing each neuron in the network with a random vm. Run the simulation with random_vm_init=True to see how the synchronization emerges over time.
../_images/Brunel_Synchronization.png

Simulation of a network with random v_m initialization. The synchronization of the neurons is not an artifact of shared initial conditions, but emerges over time.

10.4. Bonus: Power Spectrum of the Population Activity

We can get more insights into the statistics of the network activity by analysing the power spectrum of the spike trains and the population activity. The four regimes (SR, AI, SI-fast, SI-slow) are characterized by two properties: the regularity/irregularity of individual neuron’s spike trains and the stationary/oscillatory pattern of the population activity A(t). We transform the spike trains and \(A(t)\) into the frequency domain to identify regularities.

10.4.1. Question: Sampling the Population Activity

  • When analysing the population activity \(A(t)\), what is the lowest/highest frequency we are interested in?

The highest frequency \(f_{max}\) one can resolve from the time series \(A(t)\) is determined by \(\Delta t\). Even if we are not interested in very high frequencies, we should not increase \(\Delta t\) (too much) because it may affect the accuracy of the simulation.

The lowest frequency \(\Delta f\) is determined by the signal length \(T_{simulation}\). We could therefore decrease the simulation duration if we accept decreasing the resolution in the frequency domain. But there is another option: We use a “too long” simulation time \(T_{simulation}\) but then split the RateMonitor.rate signal into \(k\) chunks of duration \(T_{signal}\). We can then average the power across the \(k\) repetitions. This is what the function spike_tools.get_population_activity_power_spectrum() does - we just have to calculate the parameters first:

  • Given the values \(\Delta f = 5 Hz, \Delta t = 0.1ms, T_{init}=100ms, k=5\), compute \(T_{signal}\) and \(T_{simulation}\).
(2)\[\begin{split}\begin{array}{ccll} f_{max} = \frac{f_{sampling}}{2} = \frac{1}{2 \cdot \Delta t} \\[.2cm] N \cdot \Delta t = T_{signal} \\[.2cm] 2 \cdot f_{max} = N \cdot \Delta f \\[.2cm] T_{simulation} = k \cdot T_{signal} + T_{init} \quad , \quad k \in N \\ \end{array}\end{split}\]
\(f_{Sampling}\): sampling frequency of the signal;
\(f_{max}\): highest frequency component;
\(\Delta f\): frequency resolution in fourier domain = lowest frequency component;
\(T_{Signal}\) length of the signal;
\(\Delta t\): temporal resolution of the signal;
\(N\): number of samples (same in time- and frequency- domain)
\(T_{Simulation}\): simulation time;
\(k\): number of repetitions of the signal;
\(T_{init}\): initial part of the simulation (not used for data analysis);

10.4.2. Question: Sampling a Single Neuron Spike Train

  • The sampling of each individual neuron’s spike train is different because the signal is given as a list of timestamps (SpikeMonitor.spike_trains) and needs to be transformed into a binary vector. This is done inside the function spike_tools.get_averaged_single_neuron_power_spectrum(). Read the doc to learn how to control the sampling rate.
  • The firing rate of a single neuron can be very low and very different from one neuron to another. For that reason, we do not split the spike train into \(k\) realizations but we analyse the full spike train (\(T_{simulation}-T_{init}\)). From the simulation, we get many (\(C_E + C_I\)) spike trains and we can average across a subset of neurons. Check the doc of spike_tools.get_averaged_single_neuron_power_spectrum() to learn how to control the number of neurons of this subset.

10.4.3. Question: Single Neuron activity vs. Population Activity

We can now compute and plot the power spectrum.

%matplotlib inline
from neurodynex3.brunel_model import LIF_spiking_network
from neurodynex3.tools import plot_tools, spike_tools
import brian2 as b2

# Specify the parameters of the desired network state (e.g. SI fast)
poisson_rate = XXXX *b2.Hz
g = XXXX
CE = XXXX

# Specify the signal and simulation properties:
delta_t = XXXX * b2.ms
delta_f = XXXX * b2.Hz
T_init = XXXX * b2.ms
k = XXXX

# compute the remaining values:
f_max = XXXX
N_samples = XXXX
T_signal = XXXX
T_sim = k * T_signal + T_init

# replace the XXXX by appropriate values:

print("Start simulation. T_sim={}, T_signal={}, N_samples={}".format(T_sim, T_signal, N_samples))
b2.defaultclock.dt = delta_t
# for technical reason (solves rounding issues), we add a few extra samples:
stime = T_sim + (10 + k) * b2.defaultclock.dt
rate_monitor, spike_monitor, voltage_monitor, monitored_spike_idx = \
    LIF_spiking_network.simulate_brunel_network(
        N_Excit=CE, poisson_input_rate=poisson_rate, g=g, sim_time=stime)

plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor,
                                 spike_train_idx_list=monitored_spike_idx, t_min=0*b2.ms)
plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor,
                                 spike_train_idx_list=monitored_spike_idx, t_min=T_sim - XXXX *b2.ms)
spike_stats = spike_tools.get_spike_train_stats(spike_monitor, window_t_min= T_init)
plot_tools.plot_ISI_distribution(spike_stats, hist_nr_bins= XXXX, xlim_max_ISI= XXXX *b2.ms)

#  Power Spectrum
pop_freqs, pop_ps, average_population_rate = \
    spike_tools.get_population_activity_power_spectrum(
        rate_monitor, delta_f, k, T_init)
plot_tools.plot_population_activity_power_spectrum(pop_freqs, pop_ps, XXXX *b2.Hz, average_population_rate)
freq, mean_ps, all_ps, mean_firing_rate, all_mean_firing_freqs = \
    spike_tools.get_averaged_single_neuron_power_spectrum(
        spike_monitor, sampling_frequency=1./delta_t, window_t_min= T_init,
        window_t_max=T_sim, nr_neurons_average= XXXX )
plot_tools.plot_spike_train_power_spectrum(freq, mean_ps, all_ps, max_freq= XXXX * b2.Hz,
                                           mean_firing_freqs_per_neuron=all_mean_firing_freqs,
                                           nr_highlighted_neurons=2)
print("done")

The figures below show the type of analysis you can do with this script. The first figure shows the last 80ms of a network simulation, the second figure the power spectrum of the population activity \(A(t)\) and the third figure shows the power spectrum of single neurons (individual neurons and averaged across neurons). Note the qualitative difference between the spectral density of the population and that of the individual neurons.

../_images/Brunel_SIfast_activity.png
../_images/Brunel_SIfast_PSpop.png
../_images/Brunel_SIfast_PSsingle.png

Single neurons (red, grey) fire irregularly (I) while the population activity is synchronized (S) and oscillates. The network is in the SI regime.