In [1]:

    

### Code written for STDP tutorial @ CAMP 2016
### Adapted from Song et al 2000 (Competitive Hebbian Learning through spike-timing-dependent synaptic plasticity)
### With independent Poisson inputs (relevant figures : Fig 2)

### Author: Harsha Gurnani
### Date: June 15, 2016
#######################################


from brian2 import *
from time import time
%matplotlib inline

#### Simulation parameters ####
simtime = 150 *second
delta = 0.1*ms

defaultclock.dt = delta


    

In [2]:

    

### Neuron and synapses:
#######################
taum = 20*ms                # Membrane time constant
Vrest = -70*mV              # Resting membrane potential
Erev = 0*mV                 # Excitatory synapse - reversal potential
Irev = -70*mV               # Inhibitory synapse - reversal potential
taue = 5*ms                 # Excitatory time constants
taui = 5*ms                 # Inhibitory time constants

gmax = 0.015                #Max excitatory conductance //in unnits of leak conductance
ginh = 0.05                 #Inhibitory conductance

Vt = -54*mV                 # Spike threshold
Vr = -60*mV                 # Reset potential


    

In [3]:

    

#### How does no. of synapses/ ratio influence firing rates??
Ne = 1000                   #No. of excitatory synapses
Ni = 200                    #No. of ihibitory synapses

# How does final distribution of weights depend on presynaptic firing rates? Why???
FRe = 15*Hz                 #Firing rate (FR) for excitatory input
FRi = 10*Hz                 #FR for inhibitory input


    

In [4]:

    

### Neuron model
################

eqs = Equations('''

    dV/dt = ( (Vrest - V) + ge*(Erev - V) + gi*(Irev - V) )/taum :volt
    dge/dt = -ge/taue :1        #Conductance of exc synapse
    dgi/dt = -gi/taui :1        #Conductance of inh synapse
''')

NRN = NeuronGroup(1, model=eqs, threshold = 'V>Vt', reset = 'V=Vr', method='euler')


    

In [5]:

    

## STDP parameters for excitatory synapses
###########################################

taupre = 20*ms          # Time constant of potentiation for pre-post pairing
taupost = 20*ms         # Time constant of depression for post-pre pairing  sudo apt-get install moose moogli
Weakratio = 1.05        # Apost*taupost/(Apre*taupre)
                        # Depression:Potentiation ratio (slightly greater than 1 for stabilising network)
Apre = 0.005            # %Strengthening with pre-post pair
Apost = Apre*(taupre/taupost)*Weakratio


    

In [6]:

    

### Excitatory Synapse Model:
#############################

'''
Synapse has a synaptic weight (gsyn) and the variables x and y.
Pre: What happens when a presynaptic spike arrives at the synapse? (Depression)
Post: What happens when the postsynaptic cell spikes? (Potentiation)
'''


    

In [7]:

    

### Poisson Input at Excitatory synapses, at firing rate FRe
InpE = PoissonGroup(Ne, rates = FRe)

S_exc = Synapses(InpE, NRN, model=syneqs, on_pre=preeqs, on_post=posteqs )
S_exc.connect()                     # Connect all
S_exc.gsyn[:] = gmax*rand(Ne)       # Initialise uniformly between 0 and gmax


    

In [8]:

    

### Inhibitory synapses:
InpI = PoissonGroup(Ni, rates = FRi)

# Constant inhibitory conductance = ginh
S_inh = Synapses(InpI, NRN, model = ' ', on_pre = '''gi_post += ginh''')
S_inh.connect()


    

In [9]:

    

VMon = StateMonitor( NRN, 'V', record = 0 )                 ## Membrane potential
FR = PopulationRateMonitor(NRN)                             ## Firing rate of NRN
Weights = StateMonitor( S_exc, 'gsyn', record = [0,Ne-1])   ## Two example synaptic weights

run( simtime, report='text')


    

    




Starting simulation at t=0. s for a duration of 150. s
3.5731 s (2%) simulated in 10s, estimated 6m 50s remaining.
7.1571 s (4%) simulated in 20s, estimated 6m 39s remaining.
10.8555 s (7%) simulated in 30s, estimated 6m 25s remaining.
14.5528 s (9%) simulated in 40s, estimated 6m 12s remaining.
18.1698 s (12%) simulated in 50s, estimated 6m 3s remaining.
21.7277 s (14%) simulated in 1m 0s, estimated 5m 54s remaining.
25.333 s (16%) simulated in 1m 10s, estimated 5m 44s remaining.
28.9558 s (19%) simulated in 1m 20s, estimated 5m 34s remaining.
32.603 s (21%) simulated in 1m 30s, estimated 5m 24s remaining.
36.2044 s (24%) simulated in 1m 40s, estimated 5m 14s remaining.
39.8372 s (26%) simulated in 1m 50s, estimated 5m 4s remaining.
43.5428 s (29%) simulated in 2m 0s, estimated 4m 53s remaining.
47.1724 s (31%) simulated in 2m 10s, estimated 4m 43s remaining.
50.8482 s (33%) simulated in 2m 20s, estimated 4m 33s remaining.
54.5391 s (36%) simulated in 2m 30s, estimated 4m 23s remaining.
58.1893 s (38%) simulated in 2m 40s, estimated 4m 12s remaining.
61.8064 s (41%) simulated in 2m 50s, estimated 4m 3s remaining.
65.4838 s (43%) simulated in 3m 0s, estimated 3m 52s remaining.
68.9401 s (45%) simulated in 3m 10s, estimated 3m 43s remaining.
72.4396 s (48%) simulated in 3m 20s, estimated 3m 34s remaining.
75.9213 s (50%) simulated in 3m 30s, estimated 3m 25s remaining.
79.5084 s (53%) simulated in 3m 40s, estimated 3m 15s remaining.
83.1035 s (55%) simulated in 3m 50s, estimated 3m 5s remaining.
86.6258 s (57%) simulated in 4m 0s, estimated 2m 56s remaining.
90.2074 s (60%) simulated in 4m 10s, estimated 2m 46s remaining.
93.8395 s (62%) simulated in 4m 20s, estimated 2m 36s remaining.
97.3652 s (64%) simulated in 4m 30s, estimated 2m 26s remaining.
100.8763 s (67%) simulated in 4m 40s, estimated 2m 16s remaining.
104.4332 s (69%) simulated in 4m 50s, estimated 2m 7s remaining.
107.9733 s (71%) simulated in 5m 0s, estimated 1m 57s remaining.
111.5257 s (74%) simulated in 5m 10s, estimated 1m 47s remaining.
115.1799 s (76%) simulated in 5m 20s, estimated 1m 37s remaining.
118.7936 s (79%) simulated in 5m 30s, estimated 1m 27s remaining.
122.3798 s (81%) simulated in 5m 40s, estimated 1m 17s remaining.
125.9309 s (83%) simulated in 5m 50s, estimated 1m 7s remaining.
129.3903 s (86%) simulated in 6m 0s, estimated 57s remaining.
132.9026 s (88%) simulated in 6m 10s, estimated 48s remaining.
136.499 s (90%) simulated in 6m 20s, estimated 38s remaining.
140.0378 s (93%) simulated in 6m 30s, estimated 28s remaining.
143.6368 s (95%) simulated in 6m 40s, estimated 18s remaining.
147.2014 s (98%) simulated in 6m 50s, estimated 8s remaining.
150. s (100%) simulated in 6m 57s

In [12]:

    

figure(figsize=[10,8])
### Histogram of sinal synaptic weights
subplot(211)
hist(S_exc.gsyn[:] /gmax, 20)
ylabel('No. of synapses')
xlabel('Normalised synaptic weight')

### Initial membrane potential trace
subplot(223)
plot(VMon.t[0:3000] /ms, VMon.V[0,0:3000] /mV)
ylim([-80,-40])
ylabel('Membrane potential (mV)')
legend('Initial V')

### Final membrane potential trace
subplot(224)
plot(VMon.t[-3000:-1] /ms, VMon.V[0,-3000:-1] /mV)
ylim([-80,-40])
legend('Final V')


    

    
Out[12]:





<matplotlib.legend.Legend at 0x7ffb5e6ae810>






    

In [14]:

    

figure(figsize=[10,8])
### Evolution of Firing rate in time
subplot(211)
plot(FR.t /second, FR.smooth_rate(window='gaussian', width=100*ms)/Hz)
ylabel('Firing rate (Hz)')

### Evolution of two example synaptic weights
subplot(212)
plot(Weights.t /second, Weights.gsyn[0,:], 'b', label='Synapse 1')
plot(Weights.t /second, Weights.gsyn[1,:], 'g', label='Synapse N')
xlabel('Time (second)')
ylabel('Synaptic weight')

tight_layout()


    

In [ ]: