Source code for GraupnerBrunel2012_STDPfromCaPlasticity


#/**********************************************************************
#** This program is part of 'MOOSE', the
#** Messaging Object Oriented Simulation Environment.
#**           Copyright (C) 2003-2014 Upinder S. Bhalla. and NCBS
#** It is made available under the terms of the
#** GNU Lesser General Public License version 2.1
#** See the file COPYING.LIB for the full notice.
#**********************************************************************/

import moose
from pylab import *

[docs]def main(): """ Simulate a pseudo-STDP protocol and plot the STDP kernel that emerges from Ca plasticity of Graupner and Brunel 2012. Author: Aditya Gilra, NCBS, Bangalore, October, 2014. """ # ########################################### # Neuron models # ########################################### ## Leaky integrate and fire neuron Vrest = -65e-3 # V # resting potential Vt_base = -45e-3 # V # threshold Vreset = -55e-3 # V # in current steps, Vreset is same as pedestal R = 1e8 # Ohm tau = 10e-3 # s refrT = 2e-3 # s # ########################################### # Initialize neuron group # ########################################### ## two neurons: index 0 will be presynaptic, 1 will be postsynaptic network = moose.LIF( 'network', 2 ); moose.le( '/network' ) network.vec.Em = Vrest network.vec.thresh = Vt_base network.vec.refractoryPeriod = refrT network.vec.Rm = R network.vec.vReset = Vreset network.vec.Cm = tau/R network.vec.inject = 0. network.vec.initVm = Vrest ############################################# # Ca Plasticity parameters: synapses (not for ExcInhNetBase) ############################################# ### Cortical slice values -- Table Suppl 2 in Graupner & Brunel 2012 ### Also used in Higgins et al 2014 #tauCa = 22.6936e-3 # s # Ca decay time scale #tauSyn = 346.3615 # s # synaptic plasticity time scale ### in vitro values in Higgins et al 2014, faster plasticity #CaPre = 0.56175 # mM #CaPost = 1.2964 # mM ### in vivo values in Higgins et al 2014, slower plasticity ##CaPre = 0.33705 # mM ##CaPost = 0.74378 # mM #delayD = 4.6098e-3 # s # CaPre is added to Ca after this delay ## proxy for rise-time of NMDA #thetaD = 1.0 # mM # depression threshold for Ca #thetaP = 1.3 # mM # potentiation threshold for Ca #gammaD = 331.909 # factor for depression term #gammaP = 725.085 # factor for potentiation term #J = 5e-3 # V # delta function synapse, adds to Vm #weight = 0.43 # initial synaptic weight ## gammaP/(gammaP+gammaD) = eq weight w/o noise ## see eqn (22), noiseSD also appears ## but doesn't work here, ## weights away from 0.4 - 0.5 screw up the STDP rule!! #bistable = True # if bistable is True, use bistable potential for weights #noisy = False # use noisy weight updates given by noiseSD #noiseSD = 3.3501 # if noisy, use noiseSD (3.3501 from Higgins et al 2014) ######################################## ## DP STDP curve (Fig 2C) values -- Table Suppl 1 in Graupner & Brunel 2012 tauCa = 20e-3 # s # Ca decay time scale tauSyn = 150.0 # s # synaptic plasticity time scale CaPre = 1.0 # arb CaPost = 2.0 # arb delayD = 13.7e-3 # s # CaPre is added to Ca after this delay # proxy for rise-time of NMDA thetaD = 1.0 # mM # depression threshold for Ca thetaP = 1.3 # mM # potentiation threshold for Ca gammaD = 200.0 # factor for depression term gammaP = 321.808 # factor for potentiation term J = 5e-3 # V # delta function synapse, adds to Vm weight = 0.5 # initial synaptic weight # gammaP/(gammaP+gammaD) = eq weight w/o noise # see eqn (22), noiseSD also appears # but doesn't work here, # weights away from 0.4 - 0.5 screw up the STDP rule!! bistable = True # if bistable is True, use bistable potential for weights noisy = False # use noisy weight updates given by noiseSD noiseSD = 2.8284 # if noisy, use noiseSD (3.3501 in Higgins et al 2014) ########################################## syn = moose.GraupnerBrunel2012CaPlasticitySynHandler( '/network/syn' ) syn.numSynapses = 1 # 1 synapse # many pre-synaptic inputs can connect to a synapse # synapse onto postsynaptic neuron moose.connect( syn, 'activationOut', network.vec[1], 'activation' ) # synapse from presynaptic neuron moose.connect( network.vec[0],'spikeOut', syn.synapse[0], 'addSpike') # post-synaptic spikes also needed for STDP moose.connect( network.vec[1], 'spikeOut', syn, 'addPostSpike') syn.synapse[0].delay = 0.0 syn.synapse[0].weight = weight syn.CaInit = 0.0 syn.tauCa = tauCa syn.tauSyn = tauSyn syn.CaPre = CaPre syn.CaPost = CaPost syn.delayD = delayD syn.thetaD = thetaD syn.thetaP = thetaP syn.gammaD = gammaD syn.gammaP = gammaP syn.weightScale = J # weight ~1, weightScale ~ J # weight*weightScale is activation, # i.e. delta-fn added to postsynaptic Vm syn.weightMax = 1.0 # bounds on the weight syn.weightMin = 0. syn.noisy = noisy syn.noiseSD = noiseSD syn.bistable = bistable # ########################################### # Setting up tables # ########################################### Vms = moose.Table( '/plotVms', 2 ) moose.connect( network, 'VmOut', Vms, 'input', 'OneToOne') spikes = moose.Table( '/plotSpikes', 2 ) moose.connect( network, 'spikeOut', spikes, 'input', 'OneToOne') CaTable = moose.Table( '/plotCa', 1 ) moose.connect( CaTable, 'requestOut', syn, 'getCa') WtTable = moose.Table( '/plotWeight', 1 ) moose.connect( WtTable, 'requestOut', syn.synapse[0], 'getWeight') # ########################################### # Simulate the STDP curve with spaced pre-post spike pairs # ########################################### dt = 1e-3 # s # moose simulation moose.useClock( 0, '/network/syn', 'process' ) moose.useClock( 1, '/network', 'process' ) moose.useClock( 2, '/plotSpikes', 'process' ) moose.useClock( 3, '/plotVms', 'process' ) moose.useClock( 3, '/plotCa', 'process' ) moose.useClock( 3, '/plotWeight', 'process' ) moose.setClock( 0, dt ) moose.setClock( 1, dt ) moose.setClock( 2, dt ) moose.setClock( 3, dt ) moose.setClock( 9, dt ) moose.reinit() # function to make the aPlus and aMinus settle to equilibrium values settletime = 100e-3 # s def reset_settle(): """ Call this between every pre-post pair to reset the neurons and make them settle to rest. """ syn.synapse[0].weight = weight syn.Ca = 0.0 moose.start(settletime) # Ca gets a jump at pre-spike+delayD # So this event can occur during settletime # So set Ca and weight once more after settletime syn.synapse[0].weight = weight syn.Ca = 0.0 # function to inject a sharp current pulse to make neuron spike # immediately at a given time step def make_neuron_spike(nrnidx,I=1e-7,duration=1e-3): """ Inject a brief current pulse to make a neuron spike """ network.vec[nrnidx].inject = I moose.start(duration) network.vec[nrnidx].inject = 0. dwlist_neg = [] ddt = 2e-3 # s # since CaPlasticitySynHandler is event based # multiple pairs are needed for Ca to be registered above threshold # Values from Fig 2, last line of legend numpairs = 60 # number of spike parts per deltat t_between_pairs = 1.0 # time between each spike pair t_extent = 100e-3 # s # STDP kernel extent, # t_extent > t_between_pairs/2 inverts pre-post pairing! # dt = tpost - tpre # negative dt corresponds to post before pre print('-----------------------------------------------') for deltat in arange(t_extent,0.0,-ddt): reset_settle() for i in range(numpairs): # post neuron spike make_neuron_spike(1) moose.start(deltat) # pre neuron spike after deltat make_neuron_spike(0) moose.start(t_between_pairs) # weight changes after pre-spike+delayD # must run for at least delayD after pre-spike dw = ( syn.synapse[0].weight - weight ) / weight print(('post before pre, dt = %1.3f s, dw/w = %1.3f'%(-deltat,dw))) dwlist_neg.append(dw) print('-----------------------------------------------') # positive dt corresponds to pre before post dwlist_pos = [] for deltat in arange(ddt,t_extent+ddt,ddt): reset_settle() for i in range(numpairs): # pre neuron spike make_neuron_spike(0) moose.start(deltat) # post neuron spike after deltat make_neuron_spike(1) moose.start(t_between_pairs) dw = ( syn.synapse[0].weight - weight ) / weight print(('pre before post, dt = %1.3f s, dw/w = %1.3f'%(deltat,dw))) dwlist_pos.append(dw) print('-----------------------------------------------') print(('Each of the above pre-post pairs was repeated',\ numpairs,'times, with',t_between_pairs,'s between pairs.')) print() print('Due to event based updates, Ca decays suddenly at events:') print('pre-spike, pre-spike + delayD, and post-spike;') print('apart from the usual CaPre and CaPost jumps at') print('pre-spike + delayD and post-spike respectively.') print('Because of the event based update, multiple pre-post pairs are used.') print() print('If you reduce the t_between_pairs,') print(' you\'ll see potentiation for the LTD part without using any triplet rule!') print() print("If you turn on noise, the weights fluctuate too much,") print(" not sure if there's a bug in my noise implementation.") print('-----------------------------------------------') # ########################################### # Plot the simulated Vm-s and STDP curve # ########################################### # insert spikes so that Vm reset doesn't look weird Vmseries0 = Vms.vec[0].vector numsteps = len(Vmseries0) for t in spikes.vec[0].vector: Vmseries0[int(t/dt)-1] = 30e-3 # V Vmseries1 = Vms.vec[1].vector for t in spikes.vec[1].vector: Vmseries1[int(t/dt)-1] = 30e-3 # V timeseries = linspace(0.,1000*numsteps*dt,numsteps) # Voltage plots figure(facecolor='w') plot(timeseries,Vmseries0,color='r') # pre neuron's vm plot(timeseries,Vmseries1,color='b') # post neuron's vm xlabel('time (ms)') ylabel('Vm (V)') title("pre (r) and post (b) neurons' Vm") # Ca plots for the synapse figure(facecolor='w') plot(timeseries,CaTable.vector[:len(timeseries)],color='r') plot((timeseries[0],timeseries[-1]),(thetaP,thetaP),color='k',\ linestyle='dashed',label='pot thresh') plot((timeseries[0],timeseries[-1]),(thetaD,thetaD),color='b',\ linestyle='dashed',label='dep thresh') legend() xlabel('time (ms)') ylabel('Ca (arb)') title("Ca conc in the synapse") # Weight plots for the synapse figure(facecolor='w') plot(timeseries,WtTable.vector[:len(timeseries)],color='r') xlabel('time (ms)') ylabel('Efficacy') title("Efficacy of the synapse") # STDP curve fig = figure(facecolor='w') ax = fig.add_subplot(111) ax.plot(arange(-t_extent,0,ddt)*1000,array(dwlist_neg),'.-r') ax.plot(arange(ddt,(t_extent+ddt),ddt)*1000,array(dwlist_pos),'.-b') xmin,xmax = ax.get_xlim() ymin,ymax = ax.get_ylim() ax.set_xticks([xmin,0,xmax]) ax.set_yticks([ymin,0,ymax]) ax.plot((0,0),(ymin,ymax),linestyle='dashed',color='k') ax.plot((xmin,xmax),(0,0),linestyle='dashed',color='k') ax.set_xlabel('$t_{post}-t_{pre}$ (ms)') ax.set_ylabel('$\Delta w / w$') fig.tight_layout() #fig.subplots_adjust(hspace=0.3,wspace=0.5) # use after tight_layout() show()
if __name__ == '__main__': main()