VAbenchmarks.py

# coding: utf-8
"""
An implementation of benchmarks 1 and 2 from

    Brette et al. (2007) Journal of Computational Neuroscience 23: 349-398

The network is based on the CUBA and COBA models of Vogels & Abbott
(J. Neurosci, 2005).  The model consists of a network of excitatory and
inhibitory neurons, connected via current-based "exponential"
synapses (instantaneous rise, exponential decay).


Usage: python VAbenchmarks.py [-h] [--plot-figure] [--use-views] [--use-assembly]
                              [--use-csa] [--debug DEBUG]
                              simulator benchmark

positional arguments:
  simulator       neuron, nest, brian or another backend simulator
  benchmark       either CUBA or COBA

optional arguments:
  -h, --help      show this help message and exit
  --plot-figure   plot the simulation results to a file
  --use-views     use population views in creating the network
  --use-assembly  use assemblies in creating the network
  --use-csa       use the Connection Set Algebra to define the connectivity
  --debug DEBUG   print debugging information


Andrew Davison, UNIC, CNRS
August 2006

"""

import socket
from math import *
from pyNN.utility import get_simulator, Timer, ProgressBar, init_logging, normalized_filename
from pyNN.random import NumpyRNG, RandomDistribution


# === Configure the simulator ================================================

sim, options = get_simulator(
                    ("benchmark", "either CUBA or COBA"),
                    ("--plot-figure", "plot the simulation results to a file", {"action": "store_true"}),
                    ("--use-views", "use population views in creating the network", {"action": "store_true"}),
                    ("--use-assembly", "use assemblies in creating the network", {"action": "store_true"}),
                    ("--use-csa", "use the Connection Set Algebra to define the connectivity", {"action": "store_true"}),
                    ("--debug", "print debugging information"))

if options.use_csa:
    import csa

if options.debug:
    init_logging(None, debug=True)

timer = Timer()

# === Define parameters ========================================================

threads  = 1
rngseed  = 98765
parallel_safe = True

n        = 4000  # number of cells
r_ei     = 4.0   # number of excitatory cells:number of inhibitory cells
pconn    = 0.02  # connection probability
stim_dur = 50.   # (ms) duration of random stimulation
rate     = 100.  # (Hz) frequency of the random stimulation

dt       = 0.1   # (ms) simulation timestep
tstop    = 1000  # (ms) simulaton duration
delay    = 0.2

# Cell parameters
area     = 20000. # (µm²)
tau_m    = 20.    # (ms)
cm       = 1.     # (µF/cm²)
g_leak   = 5e-5   # (S/cm²)
if options.benchmark == "COBA":
    E_leak   = -60.  # (mV)
elif options.benchmark == "CUBA":
    E_leak   = -49.  # (mV)
v_thresh = -50.   # (mV)
v_reset  = -60.   # (mV)
t_refrac = 5.     # (ms) (clamped at v_reset)
v_mean   = -60.   # (mV) 'mean' membrane potential, for calculating CUBA weights
tau_exc  = 5.     # (ms)
tau_inh  = 10.    # (ms)

# Synapse parameters
if options.benchmark == "COBA":
    Gexc = 4.     # (nS)
    Ginh = 51.    # (nS)
elif options.benchmark == "CUBA":
    Gexc = 0.27   # (nS) #Those weights should be similar to the COBA weights
    Ginh = 4.5    # (nS) # but the delpolarising drift should be taken into account
Erev_exc = 0.     # (mV)
Erev_inh = -80.   # (mV)

### what is the synaptic delay???

# === Calculate derived parameters =============================================

area  = area*1e-8                     # convert to cm²
cm    = cm*area*1000                  # convert to nF
Rm    = 1e-6/(g_leak*area)            # membrane resistance in MΩ
assert tau_m == cm*Rm                 # just to check
n_exc = int(round((n*r_ei/(1+r_ei)))) # number of excitatory cells
n_inh = n - n_exc                     # number of inhibitory cells
if options.benchmark == "COBA":
    celltype = sim.IF_cond_exp
    w_exc    = Gexc*1e-3              # We convert conductances to uS
    w_inh    = Ginh*1e-3
elif options.benchmark == "CUBA":
    celltype = sim.IF_curr_exp
    w_exc = 1e-3*Gexc*(Erev_exc - v_mean) # (nA) weight of excitatory synapses
    w_inh = 1e-3*Ginh*(Erev_inh - v_mean) # (nA)
    assert w_exc > 0; assert w_inh < 0

# === Build the network ========================================================

extra = {'threads' : threads,
         'filename': "va_%s.xml" % options.benchmark,
         'label': 'VA'}
if options.simulator == "neuroml":
    extra["file"] = "VAbenchmarks.xml"

node_id = sim.setup(timestep=dt, min_delay=delay, max_delay=1.0, **extra)
np = sim.num_processes()

host_name = socket.gethostname()
print("Host #%d is on %s" % (node_id+1, host_name))

print("%s Initialising the simulator with %d thread(s)..." % (node_id, extra['threads']))

cell_params = {
    'tau_m'      : tau_m,    'tau_syn_E'  : tau_exc,  'tau_syn_I'  : tau_inh,
    'v_rest'     : E_leak,   'v_reset'    : v_reset,  'v_thresh'   : v_thresh,
    'cm'         : cm,       'tau_refrac' : t_refrac}

if (options.benchmark == "COBA"):
    cell_params['e_rev_E'] = Erev_exc
    cell_params['e_rev_I'] = Erev_inh

timer.start()

print("%s Creating cell populations..." % node_id)
if options.use_views:
    # create a single population of neurons, and then use population views to define
    # excitatory and inhibitory sub-populations
    all_cells = sim.Population(n_exc + n_inh, celltype(**cell_params), label="All Cells")
    exc_cells = all_cells[:n_exc]
    exc_cells.label = "Excitatory cells"
    inh_cells = all_cells[n_exc:]
    inh_cells.label = "Inhibitory cells"
else:
    # create separate populations for excitatory and inhibitory neurons
    exc_cells = sim.Population(n_exc, celltype(**cell_params), label="Excitatory_Cells")
    inh_cells = sim.Population(n_inh, celltype(**cell_params), label="Inhibitory_Cells")
    if options.use_assembly:
        # group the populations into an assembly
        all_cells = exc_cells + inh_cells

if options.benchmark == "COBA":
    ext_stim = sim.Population(20, sim.SpikeSourcePoisson(rate=rate, duration=stim_dur), label="expoisson")
    rconn = 0.01
    ext_conn = sim.FixedProbabilityConnector(rconn)
    ext_syn = sim.StaticSynapse(weight=0.1)

print("%s Initialising membrane potential to random values..." % node_id)
rng = NumpyRNG(seed=rngseed, parallel_safe=parallel_safe)
uniformDistr = RandomDistribution('uniform', low=v_reset, high=v_thresh, rng=rng)
if options.use_views:
    all_cells.initialize(v=uniformDistr)
else:
    exc_cells.initialize(v=uniformDistr)
    inh_cells.initialize(v=uniformDistr)

print("%s Connecting populations..." % node_id)
progress_bar = ProgressBar(width=20)
if options.use_csa:
    connector = sim.CSAConnector(csa.cset(csa.random(pconn)))
else:
    connector = sim.FixedProbabilityConnector(pconn, rng=rng, callback=progress_bar)
exc_syn = sim.StaticSynapse(weight=w_exc, delay=delay)
inh_syn = sim.StaticSynapse(weight=w_inh, delay=delay)

connections = {}
if options.use_views or options.use_assembly:
    connections['exc'] = sim.Projection(exc_cells, all_cells, connector, exc_syn, receptor_type='excitatory')
    connections['inh'] = sim.Projection(inh_cells, all_cells, connector, inh_syn, receptor_type='inhibitory')
    if (options.benchmark == "COBA"):
        connections['ext'] = sim.Projection(ext_stim, all_cells, ext_conn, ext_syn, receptor_type='excitatory')
else:
    connections['e2e'] = sim.Projection(exc_cells, exc_cells, connector, exc_syn, receptor_type='excitatory')
    connections['e2i'] = sim.Projection(exc_cells, inh_cells, connector, exc_syn, receptor_type='excitatory')
    connections['i2e'] = sim.Projection(inh_cells, exc_cells, connector, inh_syn, receptor_type='inhibitory')
    connections['i2i'] = sim.Projection(inh_cells, inh_cells, connector, inh_syn, receptor_type='inhibitory')
    if (options.benchmark == "COBA"):
        connections['ext2e'] = sim.Projection(ext_stim, exc_cells, ext_conn, ext_syn, receptor_type='excitatory')
        connections['ext2i'] = sim.Projection(ext_stim, inh_cells, ext_conn, ext_syn, receptor_type='excitatory')

# === Setup recording ==========================================================
print("%s Setting up recording..." % node_id)
if options.use_views or options.use_assembly:
    all_cells.record('spikes')
    exc_cells[[0, 1]].record('v')
else:
    exc_cells.record('spikes')
    inh_cells.record('spikes')
    exc_cells[0, 1].record('v')

buildCPUTime = timer.diff()

# === Save connections to file =================================================

#for prj in connections.keys():
    #connections[prj].saveConnections('Results/VAbenchmark_%s_%s_%s_np%d.conn' % (benchmark, prj, options.simulator, np))
saveCPUTime = timer.diff()

# === Run simulation ===========================================================

print("%d Running simulation..." % node_id)

sim.run(tstop)

simCPUTime = timer.diff()

E_count = exc_cells.mean_spike_count()
I_count = inh_cells.mean_spike_count()

# === Print results to file ====================================================

print("%d Writing data to file..." % node_id)

filename = normalized_filename("Results", "VAbenchmarks_%s_exc" % options.benchmark, "pkl",
                               options.simulator, np)
exc_cells.write_data(filename,
                     annotations={'script_name': __file__})
inh_cells.write_data(filename.replace("exc", "inh"),
                     annotations={'script_name': __file__})

writeCPUTime = timer.diff()

if options.use_views or options.use_assembly:
    connections = "%d e→e,i  %d i→e,i" % (connections['exc'].size(),
                                          connections['inh'].size())
else:
    connections = u"%d e→e  %d e→i  %d i→e  %d i→i" % (connections['e2e'].size(),
                                                       connections['e2i'].size(),
                                                       connections['i2e'].size(),
                                                       connections['i2i'].size())

if node_id == 0:
    print("\n--- Vogels-Abbott Network Simulation ---")
    print("Nodes                  : %d" % np)
    print("Simulation type        : %s" % options.benchmark)
    print("Number of Neurons      : %d" % n)
    print("Number of Synapses     : %s" % connections)
    print("Excitatory conductance : %g nS" % Gexc)
    print("Inhibitory conductance : %g nS" % Ginh)
    print("Excitatory rate        : %g Hz" % (E_count*1000.0/tstop,))
    print("Inhibitory rate        : %g Hz" % (I_count*1000.0/tstop,))
    print("Build time             : %g s" % buildCPUTime)
    #print("Save connections time  : %g s" % saveCPUTime)
    print("Simulation time        : %g s" % simCPUTime)
    print("Writing time           : %g s" % writeCPUTime)


# === Finished with simulator ==================================================

sim.end()