We simulate a full-scale spiking model of the local cortical microcircuit by Potjans, Tobias C., and Markus Diesmann. "The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model." Cerebral Cortex (2014): bhs358.
The model involves 4 cortical layers, each containing an axcitatory and inhibitory neuronal population:
The model also includes thalamic inputs to the cortical layers
In [1]:
import sys
sys.path.append('/opt/lib/python2.7/site-packages/')
In [2]:
import numpy as np
import pylab
In [3]:
import math
In [4]:
import nest
import nest.raster_plot
In [5]:
import logging as log
Area of network in $$mm^2$$ scales numbers of neurons
use 1 for the full-size network (77,169 neurons)
In [6]:
area = 0.1
In [7]:
preserve_K = True
In [8]:
layers = ['L23', 'L4', 'L5', 'L6']
In [9]:
full_scale_n_neurons = np.array([[20683, # layer 2/3 e
5834], # layer 2/3 i
[21915, # layer 4 e
5479], # layer 4 i
[4850, # layer 5 e
1065], # layer 5 i
[14395, # layer 6 e
2948]] # layer 6 i
)
Mean EPSP amplitude (mV) for all connections except L4e->L2/3e
In [10]:
PSP_e = .15
Mean EPSP amplitude (mv) for L4e->L2/3e connections see p. 801 of the paper, second paragraph under 'Model Parameterization', and the caption to Supplementary Fig. 7
In [11]:
PSP_e_23_4 = PSP_e * 2
Standard deviation of PSC amplitudes relative to mean PSC amplitudes
In [12]:
PSC_rel_sd = 0.1
IPSP amplitude relative to EPSP amplitude
In [13]:
g = -4.
Connection probabilities for >=1 connection between neurons in the given populations columns correspond to source populations; rows to target populations
In [14]:
# source 2/3e 2/3i 4e 4i 5e 5i 6e 6i
conn_probs = np.array([[0.1009, 0.1689, 0.0437, 0.0818, 0.0323, 0., 0.0076, 0. ], # 2/3e
[0.1346, 0.1371, 0.0316, 0.0515, 0.0755, 0., 0.0042, 0. ], # 2/3i
[0.0077, 0.0059, 0.0497, 0.135, 0.0067, 0.0003, 0.0453, 0. ], # 4e
[0.0691, 0.0029, 0.0794, 0.1597, 0.0033, 0., 0.1057, 0. ], # 4i
[0.1004, 0.0622, 0.0505, 0.0057, 0.0831, 0.3726, 0.0204, 0. ], # 5e
[0.0548, 0.0269, 0.0257, 0.0022, 0.06, 0.3158, 0.0086, 0. ], # 5i
[0.0156, 0.0066, 0.0211, 0.0166, 0.0572, 0.0197, 0.0396, 0.2252], # 6e
[0.0364, 0.001, 0.0034, 0.0005, 0.0277, 0.008, 0.0658, 0.1443]] # 6i
)
Mean dendritic delays for excitatory and inhibitory transmission (ms)
In [15]:
delays = np.array([1.5, 0.75])
Standard deviation relative to mean delays
In [16]:
delay_rel_sd = 0.5
Connection pattern used in connection calls connecting populations
In [17]:
conn_dict = {'rule': 'fixed_total_number'}
Weight distribution of connections between populations
In [18]:
weight_dict_exc = {'distribution': 'normal_clipped', 'low': 0.0}
weight_dict_inh = {'distribution': 'normal_clipped', 'high': 0.0}
Delay distribution of connections between populations
In [19]:
delay_dict = {'distribution': 'normal_clipped', 'low': 0.1}
Default synapse dictionary
In [20]:
syn_dict = {'model': 'static_synapse'}
In [21]:
neuron_model = "iaf_psc_exp"
In [22]:
Vm0_mean = -58.0
In [23]:
Vm0_std = 10.0
Neuron model parameters
In [24]:
model_params = {'tau_m': 10., # Membrane time constant (ms)
'tau_syn_ex': 0.5, # Excitatory synaptic time constant (ms)
'tau_syn_in': 0.5, # Inhibitory synaptic time constant (ms)
't_ref': 2., # Absolute refractory period (ms)
'E_L': -65., # Resting membrane potential (mV)
'V_th': -50., # Spike threshold (mV)
'C_m': 250., # Membrane capacitance (pF)
'V_reset': -65. # Reset potential (mV)
}
Rate of background Poisson input at each external input synapse (spikes/s)
In [25]:
bg_rate = 8.
DC amplitude at each external input synapse (pA) This is relevant for reproducing Potjans & Diesmann (2012) Fig. 7.
In [26]:
dc_amplitude = 0.
In-degrees for background input
In [27]:
K_bg = np.array([[1600, # 2/3e
1500], # 2/3i
[2100, # 4e
1900], # 4i
[2000, # 5e
1900], # 5i
[2900, # 6e
2100]] # 6i
)
Optional additional thalamic input (Poisson) Set n_thal to 0 to avoid this input. For producing Potjans & Diesmann (2012) Fig. 10, n_thal=902 was used. Note that the thalamic rate here reproduces the simulation results Shown in the paper, and differs from the rate given in the text.
In [28]:
n_thal = 0. # Size of thalamic population
th_start = 700. # Onset of thalamic input (ms)
th_duration = 10. # Duration of thalamic input (ms)
th_rate = 120. # Rate of thalamic neurons (spikes/s)
PSP_ext = 0.15 # Mean EPSP amplitude (mV) for external input
Connection probabilities for thalamic input
In [29]:
C_th = np.array([[0.0, # 2/3e
0.0], # 2/3i
[0.0983, # 4e
0.0619], # 4i
[0.0, # 5e
0.0], # 5i
[0.0512, # 6e
0.0196]] # 6i
)
Mean delay of thalamic input (ms)
In [30]:
delay_th = 1.5
Standard deviation relative to mean delay of thalamic input
In [31]:
delay_th_rel_sd = 0.5
In [32]:
t_sim = 1000.0 # simulated time (ms)
dt = 0.1 # simulation step (ms). ault is 0.1 ms.
allgather = True # communication protocol
In [33]:
master_seed = 123456 # Changes rng_seeds and grng_seed
n_mpi_procs = 1 # Number of MPI processes
n_threads_per_proc = 8 # Number of threads per MPI process
# Use for instance 24 for a full-scale simulation
In [34]:
n_vp = n_threads_per_proc * n_mpi_procs
In [35]:
overwrite_existing_files = True
Whether to record spikes from a fixed fraction of neurons in each population If false, a fixed number of neurons is recorded in each population. Record_fraction_neurons_spikes True with f_rec_spikes 1. records all spikes
In [36]:
record_fraction_neurons_spikes = True
if record_fraction_neurons_spikes:
frac_rec_spikes = 0.1
else:
n_rec_spikes = 100
Whether to record voltage from a fixed fraction of neurons in each population
In [37]:
record_fraction_neurons_voltage = True
if record_fraction_neurons_voltage:
frac_rec_voltage = 0.02
else:
n_rec_voltage = 20.
Whether to write any recorded cortical spikes to file
In [38]:
save_cortical_spikes = True
Whether to write any recorded membrane potentials to file
In [39]:
save_voltages = True
Whether to record thalamic spikes (only used when n_thal in Network_params.sli is nonzero)
In [40]:
record_thalamic_spikes = True
Whether to write any recorded thalamic spikes to file
In [41]:
save_thalamic_spikes = True
Name of file to which to write global IDs
In [42]:
GID_filename = 'population_GIDs.dat'
Stem for spike detector file labels
In [43]:
spike_detector_label = 'spikes'
Stem for voltmeter file labels
In [44]:
voltmeter_label = 'voltages'
Stem for thalamic spike detector file labels
In [45]:
th_spike_detector_label = 'th_spikes'
In [46]:
n_layers = None
n_pops_per_layer = None
normal_rdvs = None
PSC_e = None
PSC_e_23_4 = None
PSP_i = None
PSC_i = None
PSC_ext = None
PSC_array = None
PSC_sd = None
PSC_th_sd = None
delays_sd = None
delay_th_sd = None
n_neurons_rec_spikes = None
n_neurons_rec_voltage = None
n_neurons = None
neuron_subnet_GIDs = None
spike_detector_GIDs = None
voltmeter_GIDs = None
poisson_GIDs = None
dc_GIDs = None
th_neuron_subnet_GID = None
th_poisson_GID = None
th_spike_detector_GID = None
In [47]:
def GetLocalNodes(subnets):
if type(subnets) is not tuple:
subnets = tuple(subnets)
return nest.GetNodes(subnets, local_only=True)
In [48]:
def GetGlobalNodes(subnets):
if type(subnets) is not tuple:
subnets = tuple(subnets)
return nest.GetNodes(subnets, local_only=False)
In [49]:
def CheckParameters():
global n_layers
global n_pops_per_layer
if neuron_model != 'iaf_psc_exp':
if nest.Rank() != 0:
log.warn('Unexpected neuron type: script is tuned to "iaf_psc_exp" neurons.')
# number of layers
n_layers = len(full_scale_n_neurons)
# number of populations in each layer
n_pops_per_layer = np.shape(full_scale_n_neurons)[1]
# if np.shape(conn_probs)[0] != n_layers*n_pops_per_layer or \
# np.shape(conn_probs)[1] != n_layers*n_pops_per_layer:
# raise ValueError('conn_probs_dimensions')
# Fraction of neurons to be recorded from cannot be greater than 1
if record_fraction_neurons_spikes:
if frac_rec_spikes > 1:
raise ValueError('frac_rec_spikes')
else:
# Calculate the number of neurons from which spikes are to be recorded
if n_rec_spikes > area * min(map(min, full_scale_n_neurons)):
raise ValueError('n_rec_spikes')
# Fraction of neurons to be recorded from cannot be greater than 1
if record_fraction_neurons_voltage:
if frac_rec_voltage > 1:
raise ValueError('frac_rec_voltage')
else:
# Calculate the number of neurons from which voltages are to be recorded
if n_rec_voltage > area * min(map(min, full_scale_n_neurons)):
raise ValueError('n_rec_voltage')
In [50]:
def PrepareSimulation():
global normal_rdvs
nest.ResetKernel()
nest.SetKernelStatus({
'resolution': dt,
'total_num_virtual_procs': n_vp,
'communicate_allgather': allgather,
'overwrite_files': overwrite_existing_files,
'rng_seeds': range(master_seed, master_seed + n_vp), # local RNG seeds
'grng_seed': master_seed + n_vp # global RNG seed
})
# if run_mode == 'production':
# nest.SetKernelStatus({'data_path': output_path})
seed_offset = master_seed + n_vp
normal_rdvs = [np.random.RandomState(s) for s in range(seed_offset, seed_offset + n_vp)]
In [51]:
def DerivedParameters():
global PSC_e
global PSC_e_23_4
global PSP_i
global PSC_i
global PSC_ext
global PSC_array
global PSC_sd
global PSC_th_sd
global delays_sd
global delay_th_sd
global n_neurons_rec_spikes
global n_neurons_rec_voltage
global n_neurons
# compute numbers of neurons for the given surface area
n_neurons = np.array(map(lambda x: map(int, x), full_scale_n_neurons*area))
m = model_params
# compute PSC amplitude from PSP amplitude
# factor for transforming PSP amplitude to PSC amplitude
re = m['tau_m'] / m['tau_syn_ex']
de = m['tau_syn_ex'] - m['tau_m']
ri = m['tau_m'] / m['tau_syn_in']
di = m['tau_syn_in'] - m['tau_m']
PSC_e_over_PSP_e = (((m['C_m'])**(-1)*m['tau_m']*m['tau_syn_ex']/de*(re**(m['tau_m']/de)-re**(m['tau_syn_ex']/de)))**(-1))
PSC_i_over_PSP_i = (((m['C_m'])**(-1)*m['tau_m']*m['tau_syn_in']/di*(ri**(m['tau_m']/di)-ri**(m['tau_syn_in']/di)))**(-1))
PSC_e = PSC_e_over_PSP_e * PSP_e
PSC_e_23_4 = PSC_e_over_PSP_e * PSP_e_23_4
PSP_i = PSP_e * g
PSC_i = PSC_i_over_PSP_i * PSP_i
# PSC amplitude for all external input
PSC_ext = PSC_e_over_PSP_e * PSP_ext
# array of synaptic current amplitudes
PSC_array = np.tile(np.array([PSC_e, PSC_i]), (4,2,4,1))
PSC_array[0, 0, 1, 0] = PSC_e_23_4
# standard deviations of synaptic current amplitudes
PSC_sd = np.array([PSC_e, PSC_i]) * PSC_rel_sd
PSC_th_sd = PSC_ext * PSC_rel_sd
# standard deviations of delays
delays_sd = delays * delay_rel_sd
delay_th_sd = delay_th * delay_th_rel_sd
# numbers of neurons from which to record spikes and membrane potentials
if record_fraction_neurons_spikes:
n_neurons_rec_spikes = frac_rec_spikes*n_neurons
else:
n_neurons_rec_spikes = np.tile(n_rec_spikes, (n_layers, n_pops_per_layer, 1))
if record_fraction_neurons_voltage:
n_neurons_rec_voltage = frac_rec_voltage*n_neurons
else:
n_neurons_rec_voltage = np.tile(n_rec_voltage, (n_layers, n_pops_per_layer, 1))
In [52]:
def CreateNetworkNodes():
global neuron_subnet_GIDs
global spike_detector_GIDs
global voltmeter_GIDs
global poisson_GIDs
global dc_GIDs
global th_neuron_subnet_GID
global th_poisson_GID
global th_spike_detector_GID
# create and configure neurons
nest.SetDefaults(neuron_model, model_params)
# arrays of GIDs:
# neuron subnets
neuron_subnet_GIDs = np.tile(0, (n_layers, n_pops_per_layer, 1))
# spike detectors
spike_detector_GIDs = np.tile(0, (n_layers, n_pops_per_layer, 1))
# voltmeters
voltmeter_GIDs = np.tile(0, (n_layers, n_pops_per_layer, 1))
# Poisson generators
poisson_GIDs = np.tile(0, (n_layers, n_pops_per_layer, 1))
# DC generators
dc_GIDs = np.tile(0, (n_layers, n_pops_per_layer, 1))
for layer_index in xrange(n_layers):
nest.ChangeSubnet((0,)) # change to the root node
layer_subnet = nest.Create('subnet')
for population_index in xrange(n_pops_per_layer):
nest.ChangeSubnet(layer_subnet)
population_subnet = nest.Create('subnet')
nest.ChangeSubnet(population_subnet)
# create neurons
neuron_subnet = nest.Create('subnet')
nest.ChangeSubnet(neuron_subnet)
neuron_subnet_GIDs[layer_index][population_index] = neuron_subnet
nest.Create(neuron_model, n_neurons[layer_index][population_index])
# initialize membrane potentials
ctr = 0
for n in GetLocalNodes(neuron_subnet)[0]:
nest.SetStatus((n,), {'V_m': normal_rdvs[ctr].normal(Vm0_mean, Vm0_std)})
nest.ChangeSubnet(population_subnet)
# create and configure stimulus and recording devices
device_subnet = nest.Create('subnet')
nest.ChangeSubnet(device_subnet)
this_spike_detector = nest.Create('spike_detector')
# Set spike detector label for filenames. The GID of the spike
# detector and the process number are appended automatically.
nest.SetStatus(this_spike_detector, {
'label': spike_detector_label + '_' + str(layer_index) + '_' + str(population_index),
'to_file': save_cortical_spikes
})
spike_detector_GIDs[layer_index][population_index] = this_spike_detector
this_voltmeter = nest.Create('voltmeter')
nest.SetStatus(this_voltmeter, {
'label': voltmeter_label + '_' + str(layer_index) + '_' + str(population_index),
'to_file': save_voltages
})
voltmeter_GIDs[layer_index][population_index] = this_voltmeter
this_poisson_generator = nest.Create('poisson_generator')
this_K_bg = K_bg[layer_index][population_index]
nest.SetStatus(this_poisson_generator, {
'rate': this_K_bg * bg_rate
})
poisson_GIDs[layer_index][population_index] = this_poisson_generator
this_dc_generator = nest.Create('dc_generator')
nest.SetStatus(this_dc_generator, {
'amplitude': this_K_bg * dc_amplitude
})
dc_GIDs[layer_index][population_index] = this_dc_generator
# create and configure thalamic neurons (parrots) and their Poisson inputs
nest.ChangeSubnet((0,))
if n_thal > 0:
th_subnet = nest.Create('subnet')
nest.ChangeSubnet(th_subnet)
th_neuron_subnet_GID = nest.Create('subnet')
nest.ChangeSubnet(th_neuron_subnet_GID)
nest.Create('parrot_neuron')
nest.ChangeSubnet(th_subnet)
th_device_subnet = nest.Create('subnet')
nest.ChangeSubnet(th_device_subnet)
th_poisson_GID = nest.Create('poisson_generator')
nest.SetStatus(th_poisson_GID, {
'rate': th_rate,
'start': th_start,
'stop': th_start + th_duration
})
if record_thalamic_spikes:
th_spike_detector_GID = nest.Create('spike_detector')
nest.SetStatus(th_spike_detector_GID, {
'label': th_spike_detector_label,
'to_file': save_thalamic_spikes
})
In [53]:
def WriteGIDstoFile():
# if run_mode == 'test':
f = GID_filename
# if run_mode == 'production':
# f = output_path + '/' + GID_filename
with open(f, 'w') as f:
for n in neuron_subnet_GIDs.flatten():
GIDs = nest.GetNodes((n,))
f.write(str(min(GIDs[0])) + '\t' + str(max(GIDs[0])) + '\n')
f.close()
In [54]:
def ConnectNetworkNodes():
global neuron_subnet_GIDs
global spike_detector_GIDs
global voltmeter_GIDs
global poisson_GIDs
global dc_GIDs
global th_neuron_subnet_GID
global th_poisson_GID
global th_spike_detector_GID
# target layer
for target_layer in xrange(n_layers):
for target_pop in xrange(n_pops_per_layer):
# get neuron IDs
target_nodes = GetGlobalNodes(neuron_subnet_GIDs[target_layer][target_pop])
n_targets = n_neurons[target_layer][target_pop]
full_scale_n_targets = full_scale_n_neurons[target_layer][target_pop]
# source layer
for source_layer in xrange(n_layers):
# source population
for source_pop in xrange(n_pops_per_layer):
### local connections
# get neuron IDs
source_nodes = GetGlobalNodes(neuron_subnet_GIDs[source_layer][source_pop])
n_sources = n_neurons[source_layer][source_pop]
full_scale_n_sources = full_scale_n_neurons[source_layer][source_pop]
# get connection probability
# pick row (target) in conn_probs
r = (target_layer * n_pops_per_layer) + target_pop
# pick column (source) in conn_probs
c = (source_layer * n_pops_per_layer) + source_pop
this_conn = conn_probs[r][c]# probability for this connection
# Compute numbers of synapses assuming binomial degree
# distributions and allowing for multapses (see Potjans
# and Diesmann 2012 Cereb Cortex Eq. 1)
if preserve_K:
prod = full_scale_n_sources * full_scale_n_targets
n_syn_temp = np.log(1.-this_conn)/np.log((prod-1.)/prod)
this_n_synapses = int((n_syn_temp * n_targets) / full_scale_n_targets)
else:
prod = n_sources * n_targets
this_n_synapses = int(np.log(1.-this_conn)/np.log((prod-1.)/prod))
if this_n_synapses > 0:
mean_weight = PSC_array[target_layer][target_pop][source_layer][source_pop]
# Create label for target and source populations
conn_label = 'layers' + str(target_layer) + '_' + 'populations' + str(target_pop) + '-' + \
'layers' + str(source_layer) + '_' + 'populations' + str(source_pop)
# fill the weight dictionary for Connect and insert it into the synapse dictionary
if mean_weight > 0:
weight_dict_exc['mu'] = mean_weight
weight_dict_exc['sigma'] = np.abs(PSC_sd[source_pop])
syn_dict['weight'] = weight_dict_exc
else:
weight_dict_inh['mu'] = mean_weight
weight_dict_inh['sigma'] = np.abs(PSC_sd[source_pop])
syn_dict['weight'] = weight_dict_inh
# fill the delay dictionary for Connect and insert it into the synapse dictionary
delay_dict['mu'] = delays[source_pop]
delay_dict['sigma'] = np.abs(delays_sd[source_pop])
syn_dict['delay'] = delay_dict
# fill the connectivity dictionary with the number of synapses to be used
conn_dict['N'] = this_n_synapses
conn_dict['rule'] = 'fixed_total_number'
# Connect the populations
nest.Connect(source_nodes[0], target_nodes[0], conn_dict, syn_dict)
if n_thal > 0:
# connections from thalamic neurons
source_nodes = GetGlobalNodes(th_neuron_subnet_GID)
this_conn = C_th[target_layer][target_pop]
# Compute numbers of synapses assuming binomial degree
# distributions and allowing for multapses (see Potjans and
# Diesmann 2012 Cereb Cortex Eq. 1)
if preserve_K:
prod = n_thal * full_scale_n_targets
n_syn_temp = np.log(1.-this_conn)/np.log((prod-1.)/prod)
this_n_synapses = int(full_scale_n_targets / (n_syn_temp * n_targets))
else:
prod = n_thal * n_targets
this_n_synapses = int(np.log(1.-this_conn)/np.log((prod-1.)/prod))
if this_n_synapses > 0:
# create label for current target population
th_conn_label = layers[target_layer] + '_' + populations[target_pop]
# fill the weight dictionary for Connect
weight_dict_exc['mu'] = PSC_ext
weight_dict_exc['sigma'] = np.abs(PSC_th_sd)
# insert the weight dictionary into the synapse dictionary
syn_dict['weight'] = weight_dict_exc
# fill the delay dictionary for Connect
delay_dict['mu'] = delay_th
delay_dict['sigma'] = np.abs(delay_th_sd)
# insert the delay dictionary into the synapse dictionary
syn_dict['delay'] = delay_dict
conn_dict['N'] = this_n_synapses
conn_dict['rule'] = 'fixed_total_number'
nest.Connect(source_nodes, target_nodes, conn_dict, syn_dict)
# Connect devices
# record from a continuous range of IDs
# (appropriate for networks without topology)
# print tuple(spike_detector_GIDs[target_layer][target_pop])
# print target_nodes[:int(n_neurons_rec_spikes[target_layer][target_pop])][0]
nest.Connect(target_nodes[:int(n_neurons_rec_spikes[target_layer][target_pop])][0],
tuple(spike_detector_GIDs[target_layer][target_pop]),
'all_to_all')
nest.Connect(tuple(voltmeter_GIDs[target_layer][target_pop]),
tuple(target_nodes[:int(n_neurons_rec_voltage[target_layer][target_pop])])[0],
'all_to_all')
nest.Connect(tuple(poisson_GIDs[target_layer][target_pop]),
target_nodes[0],
'all_to_all',
{'weight': PSC_ext, 'delay': delays[0]})
nest.Connect(tuple(dc_GIDs[target_layer][target_pop]),
target_nodes[0],
'all_to_all')
if n_thal > 0:
# Connect thalamic poisson_generator to thalamic neurons (parrots)
nest.Connect(th_poisson_GID, GetGlobalNodes(th_neuron_subnet_GID))
if record_thalamic_spikes and n_thal > 0:
# Connect thalamic neurons to spike detector
nest.Connect(GetGlobalNodes(th_neuron_subnet_GID), th_spike_detector_GID)
In [55]:
CheckParameters()
PrepareSimulation()
DerivedParameters()
CreateNetworkNodes()
WriteGIDstoFile()
ConnectNetworkNodes()
nest.Simulate(t_sim)
In [56]:
import glob
import matplotlib.pyplot as plt
import os
In [57]:
datapath = '.'
In [58]:
T = t_sim
T_start = 200. # starting point of analysis (to avoid transients)
load GIDs dumped into file
In [59]:
gidfile = open(os.path.join(datapath , 'population_GIDs.dat'), 'r')
gids = []
for l in gidfile:
a = l.split()
gids.append([int(a[0]),int(a[1])])
print 'Global IDs:'
print gids
print
calculate number of populations from the GIDs
In [60]:
num_pops = len(gids)
print 'Number of populations:'
print num_pops
print
first GID in each population
In [61]:
raw_first_gids = [gids[i][0] for i in np.arange(len(gids))]
population sizes
In [62]:
pop_sizes = [gids[i][1]-gids[i][0]+1 for i in np.arange(len(gids))]
numbers of neurons for which spikes were recorded
In [63]:
if record_fraction_neurons_spikes == True:
rec_sizes = [int(pop_sizes[i]*frac_rec_spikes) for i in xrange(len(pop_sizes))]
else:
rec_sizes = [n_rec_spikes]*len(pop_sizes)
first GID of each population once device GIDs are dropped
In [64]:
first_gids=[int(1 + np.sum(pop_sizes[:i])) for i in np.arange(len(pop_sizes))]
last GID of each population once device GIDs are dropped
In [65]:
last_gids = [int(np.sum(pop_sizes[:i+1])) for i in np.arange(len(pop_sizes))]
convert lists to a nicer format, i.e. [[2/3e, 2/3i], []....]
In [66]:
Pop_sizes =[pop_sizes[i:i+2] for i in xrange(0,len(pop_sizes),2)]
print 'Population sizes:'
print Pop_sizes
print
In [67]:
Raw_first_gids =[raw_first_gids[i:i+2] for i in xrange(0,len(raw_first_gids),2)]
First_gids = [first_gids[i:i+2] for i in xrange(0,len(first_gids),2)]
Last_gids = [last_gids[i:i+2] for i in xrange(0,len(last_gids),2)]
total number of neurons in the simulation
In [68]:
num_neurons = last_gids[len(last_gids)-1]
print 'Total number of neurons:'
print num_neurons
print
load spikes from gdf files, correct GIDs and merge them in population files, and store spike trains
In [69]:
# will contain neuron id resolved spike trains
neuron_spikes = [[] for i in np.arange(num_neurons+1)]
# container for population-resolved spike data
spike_data= [[[],[]],[[],[]],[[],[]],[[],[]],[[],[]],[[],[]],[[],[]],[[],[]]]
counter = 0
for layer in ['0','1','2','3']:
for population in ['0','1']:
output = os.path.join(datapath, 'population_spikes-{}-{}.gdf'.format(layer, population))
file_pattern = os.path.join(datapath, 'spikes_{}_{}*'.format(layer, population))
files = glob.glob(file_pattern)
print 'Merge '+str(len(files))+' spike files from L'+layer+'P'+population
if files:
merged_file = open(output,'w')
for f in files:
data = open(f,'r')
for l in data :
a = l.split()
a[0] = int(a[0])
a[1] = float(a[1])
raw_first_gid = Raw_first_gids[int(layer)][int(population)]
first_gid = First_gids[int(layer)][int(population)]
a[0] = a[0] - raw_first_gid + first_gid
if(a[1] > T_start): # discard data in the start-up phase
spike_data[counter][0].append(num_neurons-a[0])
spike_data[counter][1].append(a[1]-T_start)
neuron_spikes[a[0]].append(a[1]-T_start)
converted_line = str(a[0]) + '\t' + str(a[1]) +'\n'
merged_file.write(converted_line)
data.close()
merged_file.close()
counter +=1
Set color palette for matplotlib
In [70]:
import matplotlib.cm as cm
tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),
(44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),
(148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),
(227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),
(188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]
for i in range(len(tableau20)):
r, g, b = tableau20[i]
tableau20[i] = (r / 255., g / 255., b / 255.)
clrs = tableau20
plt.ion()
raster plot of spikes
In [71]:
plt.figure(1)
counter = 1
for j in np.arange(num_pops):
for i in np.arange(first_gids[j],first_gids[j]+rec_sizes[j]):
if len(neuron_spikes[i]):
plt.plot(neuron_spikes[i],np.ones_like(neuron_spikes[i])+sum(rec_sizes)-counter,'k o',ms=1, mfc=clrs[j],mec=clrs[j])
counter+=1
plt.xlim(0,T-T_start)
plt.ylim(0,sum(rec_sizes))
plt.xlabel(r'time (ms)')
plt.ylabel(r'neuron id')
plt.show()
Calculate firing rates
In [72]:
rates = []
temp = 0
for i in np.arange(num_pops):
for j in np.arange(first_gids[i], last_gids[i]):
temp += len(neuron_spikes[j])
rates.append(temp/(rec_sizes[i]*(T-T_start))*1e3)
temp = 0
Plot firing rates in different layers as a bar chart
In [73]:
plt.figure(2)
ticks= np.arange(num_pops)
plt.bar(ticks, rates, width=0.9, color=clrs)
xticklabels = ['L2/3e','L2/3i','L4e','L4i','L5e','L5i','L6e','L6i']
plt.setp(plt.gca(), xticks=ticks+0.5, xticklabels=xticklabels)
plt.xlabel(r'subpopulation')
plt.ylabel(r'firing rate (spikes/s)')
plt.show()
In [73]: