Source code for neurodynex3.brunel_model.LIF_spiking_network

Implementation of the Brunel 2000 network:
sparsely connected network of identical LIF neurons (Model A).

# This file is part of the exercise code repository accompanying
# the book: Neuronal Dynamics (see
# located at

# This free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License 2.0 as published by the
# Free Software Foundation. You should have received a copy of the
# GNU General Public License along with the repository. If not,
# see

# Should you reuse and publish the code for your own purposes,
# please cite the book or point to the webpage

# Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski.
# Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition.
# Cambridge University Press, 2014.

import brian2 as b2
from brian2 import NeuronGroup, Synapses, PoissonInput
from brian2.monitors import StateMonitor, SpikeMonitor, PopulationRateMonitor
from random import sample
from import plot_tools
from numpy import random
import matplotlib.pyplot as plt

# Default parameters of a single LIF neuron
V_REST = 0. * b2.mV
V_RESET = +10. * b2.mV
# Default parameters of the network
SYNAPTIC_WEIGHT_W0 = 0.1 * b2.mV
# note: w_ee = w_ei = w0 and w_ie=w_ii = -g*w0

b2.defaultclock.dt = 0.05 *

[docs]def simulate_brunel_network( N_Excit=5000, N_Inhib=None, N_extern=N_POISSON_INPUT, connection_probability=CONNECTION_PROBABILITY_EPSILON, w0=SYNAPTIC_WEIGHT_W0, g=RELATIVE_INHIBITORY_STRENGTH_G, synaptic_delay=SYNAPTIC_DELAY, poisson_input_rate=POISSON_INPUT_RATE, w_external=None, v_rest=V_REST, v_reset=V_RESET, firing_threshold=FIRING_THRESHOLD, membrane_time_scale=MEMBRANE_TIME_SCALE, abs_refractory_period=ABSOLUTE_REFRACTORY_PERIOD, monitored_subset_size=100, random_vm_init=False, sim_time=100.* """ Fully parametrized implementation of a sparsely connected network of LIF neurons (Brunel 2000) Args: N_Excit (int): Size of the excitatory popluation N_Inhib (int): optional. Size of the inhibitory population. If not set (=None), N_Inhib is set to N_excit/4. N_extern (int): optional. Number of presynaptic excitatory poisson neurons. Note: if set to a value, this number does NOT depend on N_Excit and NOT depend on connection_probability (this is different from the book and paper. Only if N_extern is set to 'None', then N_extern is computed as N_Excit*connection_probability. connection_probability (float): probability to connect to any of the (N_Excit+N_Inhib) neurons CE = connection_probability*N_Excit CI = connection_probability*N_Inhib Cexternal = N_extern w0 (float): Synaptic strength J g (float): relative importance of inhibition. J_exc = w0. J_inhib = -g*w0 synaptic_delay (Quantity): Delay between presynaptic spike and postsynaptic increase of v_m poisson_input_rate (Quantity): Poisson rate of the external population w_external (float): optional. Synaptic weight of the excitatory external poisson neurons onto all neurons in the network. Default is None, in that case w_external is set to w0, which is the standard value in the book and in the paper Brunel2000. The purpose of this parameter is to see the effect of external input in the absence of network feedback(setting w0 to 0mV and w_external>0). v_rest (Quantity): Resting potential v_reset (Quantity): Reset potential firing_threshold (Quantity): Spike threshold membrane_time_scale (Quantity): tau_m abs_refractory_period (Quantity): absolute refractory period, tau_ref monitored_subset_size (int): nr of neurons for which a VoltageMonitor is recording Vm random_vm_init (bool): if true, the membrane voltage of each neuron is initialized with a random value drawn from Uniform(v_rest, firing_threshold) sim_time (Quantity): Simulation time Returns: (rate_monitor, spike_monitor, voltage_monitor, idx_monitored_neurons) PopulationRateMonitor: Rate Monitor SpikeMonitor: SpikeMonitor for ALL (N_Excit+N_Inhib) neurons StateMonitor: membrane voltage for a selected subset of neurons list: index of monitored neurons. length = monitored_subset_size """ if N_Inhib is None: N_Inhib = int(N_Excit/4) if N_extern is None: N_extern = int(N_Excit*connection_probability) if w_external is None: w_external = w0 J_excit = w0 J_inhib = -g*w0 lif_dynamics = """ dv/dt = -(v-v_rest) / membrane_time_scale : volt (unless refractory)""" network = NeuronGroup( N_Excit+N_Inhib, model=lif_dynamics, threshold="v>firing_threshold", reset="v=v_reset", refractory=abs_refractory_period, method="linear") if random_vm_init: network.v = random.uniform(v_rest/b2.mV, high=firing_threshold/b2.mV, size=(N_Excit+N_Inhib))*b2.mV else: network.v = v_rest excitatory_population = network[:N_Excit] inhibitory_population = network[N_Excit:] exc_synapses = Synapses(excitatory_population, target=network, on_pre="v += J_excit", delay=synaptic_delay) exc_synapses.connect(p=connection_probability) inhib_synapses = Synapses(inhibitory_population, target=network, on_pre="v += J_inhib", delay=synaptic_delay) inhib_synapses.connect(p=connection_probability) external_poisson_input = PoissonInput(target=network, target_var="v", N=N_extern, rate=poisson_input_rate, weight=w_external) # collect data of a subset of neurons: monitored_subset_size = min(monitored_subset_size, (N_Excit+N_Inhib)) idx_monitored_neurons = sample(range(N_Excit+N_Inhib), monitored_subset_size) rate_monitor = PopulationRateMonitor(network) # record= some_list is not supported? :-( spike_monitor = SpikeMonitor(network, record=idx_monitored_neurons) voltage_monitor = StateMonitor(network, "v", record=idx_monitored_neurons) return rate_monitor, spike_monitor, voltage_monitor, idx_monitored_neurons
[docs]def getting_started(): """ A simple example to get started """ rate_monitor, spike_monitor, voltage_monitor, monitored_spike_idx = simulate_brunel_network( N_Excit=2000, sim_time=800. * plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min=0.*, N_highlighted_spiketrains=3, avg_window_width=1. *
def _demo_emergence_of_oscillation(): poisson_rate = 18 * b2.Hz g = 2.5 rate_monitor, spike_monitor, voltage_monitor, monitored_spike_idx = \ simulate_brunel_network(N_Excit=6000, random_vm_init=True, poisson_input_rate=poisson_rate, g=g, sim_time=300. *, monitored_subset_size=50) plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min=0* plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_max=50* plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min=250* def _some_example_calls_and_tests(): from import spike_tools poisson_rate = 35*b2.Hz g = 4 CE = 5000 delta_t = 0.1 * delta_f = 5. * b2.Hz T_init = 100 * k = 9 f_max = 1./(2. * delta_t) N_samples = 2. * f_max / delta_f T_signal = N_samples * delta_t T_sim = k * T_signal + T_init print("Start simulation. T_sim={}, T_signal={}. N_samples={}".format(T_sim, T_signal, N_samples)) b2.defaultclock.dt = delta_t stime = T_sim + (10 + k) * b2.defaultclock.dt # add a few extra samples (solves rounding issues) rate_monitor, spike_monitor, voltage_monitor, monitored_spike_idx = \ 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* plot_tools.plot_network_activity(rate_monitor, spike_monitor, voltage_monitor, spike_train_idx_list=monitored_spike_idx, t_min=T_sim - 80* spike_stats = spike_tools.get_spike_train_stats(spike_monitor, window_t_min=150.* plot_tools.plot_ISI_distribution(spike_stats, hist_nr_bins=77, xlim_max_ISI=100* # # Power Spectrum pop_freqs, pop_ps, average_population_rate = \ spike_tools.get_population_activity_power_spectrum( rate_monitor, delta_f, k, T_init, subtract_mean_activity=True) plot_tools.plot_population_activity_power_spectrum(pop_freqs, pop_ps, 1000*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=100.*, window_t_max=T_sim, subtract_mean=False, nr_neurons_average=200) print("plot_spike_train_power_spectrum") plot_tools.plot_spike_train_power_spectrum(freq, mean_ps, all_ps, 1000 * b2.Hz, mean_firing_freqs_per_neuron=all_mean_firing_freqs, nr_highlighted_neurons=2) print("done") if __name__ == "__main__": getting_started()