Introduction

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:

  1. Layer 2/3
  • Layer 4
  • Layer 5
  • Layer 6

The model also includes thalamic inputs to the cortical layers

Setup

  • Ensure that PYTHONPATH includes the nest installation directory
  • Import python numeric and plotting libraries
  • Import the python math library
  • Import nest and its required modules

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

Network parameters

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

Cortical layer population sizes

Whether to use full-scale in-degrees when downscaling the number of neurons When preserve_K is false, the full-scale connection probabilities are used.

Note that this produces different dynamics compared to the original model.


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
                                )

Connection properties

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'}

Single-neuron parameters

Neuron model. For PSP-to-PSC conversion to be correct, synapses should be current-based with an exponential time course mean of initial membrane potential (mV) std of initial membrane potential (mV)


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)
               }

Stimulus parameters

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

Simulation parameters


In [32]:
t_sim = 1000.0      # simulated time (ms)
dt = 0.1            # simulation step (ms). ault is 0.1 ms.
allgather = True    # communication protocol

Master seed for random number generators

Actual seeds will be master_seed ... master_seed + 2*n_vp ==>> different master seeds must be spaced by at least 2*n_vp + 1

See Gewaltig et al. (2012) for details


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

Number of virtual processes

This should be an integer multiple of the number of MPI processes. See Morrison et al. (2005) Neural Comput


In [34]:
n_vp = n_threads_per_proc * n_mpi_procs

Recording parameters


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'

Global variables


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

Utility functions

  1. GetLocalNodes: Fetch nodes in the local machine associated with the given subnet
  • GetGlobalNodes: Fetch all nodes in all machines associated with the given subnet
  • CheckParameters: Verify whether the parameters supplied above are sane

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')

Prepare Simulation

  1. Reset the kernel
  • Set default kernel status according to the parameters supplied above
  • Initialize thread-local random generator seeds

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)]

Calculate network and model parameters


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))

Create nodes


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
      })

Write global IDs of nodes to a file


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()

Connect the nodes


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)

Run the simulation


In [55]:
CheckParameters()

PrepareSimulation()

DerivedParameters()

CreateNetworkNodes()

WriteGIDstoFile()

ConnectNetworkNodes()

nest.Simulate(t_sim)

Analyze output


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


Global IDs:
[[4, 2071], [2079, 2661], [2670, 4860], [4868, 5414], [5423, 5907], [5915, 6020], [6029, 7467], [7475, 7768]]

calculate number of populations from the GIDs


In [60]:
num_pops = len(gids)
print 'Number of populations:'
print num_pops
print


Number of populations:
8

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


Population sizes:
[[2068, 583], [2191, 547], [485, 106], [1439, 294]]


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


Total number of neurons:
7713

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


Merge 8 spike files from L0P0
Merge 8 spike files from L0P1
Merge 8 spike files from L1P0
Merge 8 spike files from L1P1
Merge 8 spike files from L2P0
Merge 8 spike files from L2P1
Merge 8 spike files from L3P0
Merge 8 spike files from L3P1

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]: