VogelsAbbott2005 Network

sim.py

from __future__ import division
import sys
import time
from brian2 import *
#import matplotlib.pyplot as plt

delays = False
method = 'euler'

# RNG seed
seed = 2
numpy.random.seed(seed)
start_time = time.time()

prefs.codegen.target = 'numpy'

# Dynamics
eqs = Equations('''
dv/dt = (g_L*(V_L-v)+ge*(V_E-v)+gi*(V_I-v))/C : volt (unless refractory)
dge/dt = -ge*(1./tau_e) : siemens
dgi/dt = -gi*(1./tau_i) : siemens
''')

class ProgressBar(object):
    def __init__(self, toolbar_width):
        self.toolbar_width = toolbar_width
        self.ticks = 0

    def __call__(self, elapsed, complete, start, duration):
        if complete == 0.0:
            # setup toolbar
            sys.stdout.write("[%s]" % (" " * self.toolbar_width))
            sys.stdout.flush()
            sys.stdout.write("\b" * (self.toolbar_width + 1)) # return to start of line, after '['
        else:
            ticks_needed = int(round(complete * 40))
            if self.ticks < ticks_needed:
                sys.stdout.write("-" * (ticks_needed-self.ticks))
                sys.stdout.flush()
                self.ticks = ticks_needed
        if complete == 1.0:
            sys.stdout.write("\n")


class Neurons():
    def __init__(self):
        pass

    def create_group(self):
        params = self.params
        group = NeuronGroup(params['N'], model=eqs,
                    threshold='v > V_thresh',
                    reset='v = V_reset',
                    namespace=self.params,
                    refractory=5*ms,
                    method=method)

        return group


class FastSpiking(Neurons):
    def __init__(self):
        V_L = -60 * mV
        V_thresh = -50 * mV
        tau_m = 20 * ms 
        R = 100*Mohm

        self.params = {
            'g_L'      : 1/R,
            'C'        : tau_m/R,
            # Cell parameters
            'N'        : 2000,
            'V_E'      : 0 * mV,
            'V_I'      : -80 * mV,
            'V_L'      : V_L,
            'tau_m'    : tau_m,

            # Reset/threshold
            'V_thresh' : V_thresh,
            'V_reset'  : V_L,

            # Conductance time constants and amplitudes
            'tau_e'    : 5 * ms,
            'tau_i'    : 10 * ms,
            'g_bar_e'  : 3.2*nS,
            'g_bar_i'  : 52*nS
        }
        self.group = self.create_group()


class Pyramidal(Neurons):
    def __init__(self):
        V_L = -60 * mV
        V_thresh = -50 * mV
        tau_m = 20 * ms 
        R = 100*Mohm

        self.params = {
            'g_L'      : 1/R,
            'C'        : tau_m/R,
            # Cell parameters
            'N'        : 8000,
            'V_E'      : 0 * mV,
            'V_I'      : -80 * mV,
            'V_L'      : V_L,
            'tau_m'    : tau_m,

            # Reset/threshold
            'V_thresh' : V_thresh,
            'V_reset'  : V_L,

            # Conductance time constants and amplitudes
            'tau_e'    : 5 * ms,
            'tau_i'    : 10 * ms,
            'g_bar_e'  : 3.2*nS, 
            'g_bar_i'  : 52*nS  
        }
        self.group = self.create_group()
        

# Create neuron groups
Pyr = Pyramidal()
FS = FastSpiking()
E = Pyr.group
I = FS.group

# Create external Poisson units
rates = {'e_rate': 2000*Hz, 'i_rate': 2000*Hz}
Ext_E = NeuronGroup(Pyr.params['N'], 
            model='rate = e_rate : Hz', 
            threshold='rand()<(rate*dt) and (t<50*ms)', 
            namespace=rates)
Ext_I = NeuronGroup(FS.params['N'], 
            model='rate = i_rate : Hz', 
            threshold='rand()<(rate*dt) and (t<50*ms)', 
            namespace=rates)

# Initial values 
E.v = Pyr.params['V_L']
I.v = FS.params['V_L']

# Initialize synapses
Cee = Synapses(E, E, on_pre='ge += g_bar_e', namespace=Pyr.params, method=method)
Cii = Synapses(I, I, on_pre='gi += g_bar_i', namespace=FS.params, method=method)
Cie = Synapses(I, E, on_pre='gi += g_bar_i', namespace=Pyr.params, method=method)
Cei = Synapses(E, I, on_pre='ge += g_bar_e', namespace=FS.params, method=method)
Cext_e = Synapses(Ext_E, E, on_pre='ge += g_bar_e', namespace=Pyr.params)
Cext_i = Synapses(Ext_I, I, on_pre='ge += g_bar_e', namespace=FS.params)

# Connect groups
Cee.connect('i != j', p=0.02)
Cii.connect('i != j', p=0.02)
Cie.connect(True, p=0.02)
Cei.connect(True, p=0.02)
Cext_e.connect('i==j')
Cext_i.connect('i==j')

# Set delays
if delays:
    Cee.delay = 2*rand(len(Cee)) * ms
    Cii.delay = 2*rand(len(Cii)) * ms
    Cie.delay = 2*rand(len(Cie)) * ms
    Cei.delay = 2*rand(len(Cei)) * ms

# Monitor spiking
M_E = SpikeMonitor(E)
M_I = SpikeMonitor(I)
M = PopulationRateMonitor(E)

#trace_E = StateMonitor(E, ['v', 'ge', 'gi'], record=[0,1,2])
#trace_I = StateMonitor(I, ['v', 'ge', 'gi'], record=[0,1,2])

# Simulation and printout
print "Network construction time:", time.time() - start_time, "seconds"
print "Simulation running..."
start_time = time.time()
run(1*second, 
    report=ProgressBar(toolbar_width=50),
    report_period=0.5*second
)
duration = time.time() - start_time
print "Simulation time:", duration, "seconds"
print "Method:", method
print "Delays:", delays

# Average spike count across all neurons in second half of simulation
r0_E = M_E.i[M_E.t > 0.5*second].size/Pyr.params['N']/0.5
r0_I = M_I.i[M_I.t > 0.5*second].size/FS.params['N']/0.5
print r0_E, "Hz"
print r0_I, "Hz"

# Plots
#plt.subplot(611)
#plt.plot(M_E.t/ms, M_E.i, '.')
#
#plt.subplot(612)
#plt.plot(M_I.t/ms, M_I.i, '.')
#
#plt.subplot(613)
#plt.plot(trace_E.t/ms, trace_E[0].v/mV)
#
#plt.subplot(614)
#plt.plot(trace_E.t/ms, trace_E[0].ge/nS)
#
#plt.subplot(615)
#plt.plot(trace_E.t/ms, trace_E[0].gi/nS)
#
#plt.subplot(616)
#plt.plot(M.t/ms, M.rate/Hz)