######################################################
# Recreate the STDP window for pairings at 20Hz, and for both distal and proximal locations
######################################################
from __future__ import division
from brian2 import *
import numpy as np
import matplotlib.pylab as plt
import json
import os, sys
mod_path = os.path.abspath(os.path.join('..','0. Model'))
sys.path.append(mod_path)
from oo_Parameters import *
from oo_equations_AMPAplast import *
from oo_initScripts import set_init_nrn, set_init_syn
from MakeNeuron_AMPAplast import *
from MorphologyData import *
for paramS in [0,1]:
start_scope()
######################################################
## Load Morpho
######################################################
# morph = '../0. Model/Branco2010_Morpho.swc'
# morph_data = BrancoData
#
# distal_compartments_nonmda = distal_compartments_Branco_nonmda
# distal_compartments_eff = distal_compartments_Branco_eff
# proximal_compartments = proximal_compartments_Branco
morph = '../0. Model/Acker2008.swc'
morph_data = AckerData
distal_compartments_nonmda = distal_compartments_Acker_nonmda
distal_compartments_eff = distal_compartments_Acker_eff
proximal_compartments = proximal_compartments_Acker
#####################################################
# Sim parameters
#####################################################
Theta_low = morph_data['thetalow']*mV
if paramS == 0:
d_compartm = proximal_compartments
nrIn = len(d_compartm)
str_var = 'prox'
elif paramS ==1:
d_compartm = distal_compartments_eff
nrIn = len(d_compartm)
str_var = 'disteff'
elif paramS ==2:
d_compartm = distal_compartments_nonmda
nrIn = len(d_compartm)
str_var = 'distnonmda'
print('***')
if morph[12:-8] == 'Acker':
print('-- L5 '+str_var+'--')
else:
print('-- L2/3 '+str_var+'--')
pairing_rate = 20. # pairing rate
dt_array = np.array([1.,2.5,5.,7.5,10.,12.5,15.,17.5,20.]) # interspike intervals
init_weight = 0.5 # initial weight
reps = 5 # number of pairings
#####################################################
# Input neurons
#####################################################
V_rest = -70.*mV # resting potential
tau_in = 8.*ms #membrane timeconstant
V_thresh = -45.*mV # spike threshold
C = 200.*pF # membrane capacitance
# Equations input neuron
eqs_in = '''
dv/dt = (V_rest-v)/tau_in + Idrive/C: volt
Idrive : amp
ds_trace/dt = -s_trace/taux :1
'''
#####################################################
# create neuron objects
#####################################################
# IandF input neurons
N_input = NeuronGroup(2*nrIn, eqs_in, threshold='v>V_thresh',
reset='v=V_rest;s_trace+=x_reset*(taux/ms)', method='linear')#
# morphological neurons
test_model = BRIANModel(morph)
neuron = test_model.makeNeuron_Ca(morph_data)
neuron.run_regularly('Mgblock = 1./(1.+ exp(-0.062*vu2)/3.57)',dt=defaultclock.dt)
neuron2 = test_model.makeNeuron_Ca(morph_data)
neuron2.run_regularly('Mgblock = 1./(1.+ exp(-0.062*vu2)/3.57)',dt=defaultclock.dt)
print('Neurons created...')
#####################################################
# create Synapses
#####################################################
Syn_1 = Synapses(N_input,neuron,
model= eq_1_plastAMPA,
on_pre = eq_2_plastAMPA,
method='heun'
)
Syn_2 = Synapses(N_input,neuron2,
model= eq_1_plastAMPA,
on_pre = eq_2_plastAMPA,
method='heun'
)
for jj in range(nrIn):
Syn_1.connect(i=jj,j=neuron[d_compartm[jj]:d_compartm[jj]+1])
Syn_2.connect(i=nrIn+jj,j=neuron2[d_compartm[jj]:d_compartm[jj]+1])
print('Synapses created...')
for ttt in range(nrIn):
print('Start compartment '+str(ttt+1)+' of '+ str(nrIn))
#####################################################
# Set Initial Neuron Parameters
#####################################################
set_init_syn(Syn_1,init_weight)
set_init_syn(Syn_2,init_weight)
N_input.v = V_rest
N_input.s_trace = 0.
#####################################################
# Run
#####################################################
# initialize matrices to store data
weight_change1 = np.zeros(shape(dt_array))
weight_change2 = np.zeros(shape(dt_array))
print('Start simulation...')
# loop over various interspike intervals
for jj in range(size(dt_array)):
dt_var = dt_array[jj]
print('-> dt '+str(dt_var))
pair_interval = 1000./pairing_rate-3.-dt_var
set_init_syn(Syn_1,init_weight)
set_init_syn(Syn_2,init_weight)
# Initial values
set_init_nrn(neuron,Theta_low)
set_init_nrn(neuron2,Theta_low)
N_input.v = V_rest
N_input.s_trace = 0.
run(100*ms)
# loop over Pairings
for ii in range(reps):
neuron.I = 0.*pA
neuron2.I = 0.*pA
N_input.Idrive = 0.*mA
###### 1st SPIKE
neuron2.main.I = 1000.*pA
N_input.Idrive[ttt] = 2000.*pA
if dt_var >= 3:
run(3*ms)
neuron2.I = 0.*pA
N_input.Idrive = 0.*mA
run(dt_var*ms-3*ms)
neuron.main.I = 1000.*pA
N_input.Idrive[nrIn+ttt] = 2000.*pA
run(3*ms)
neuron.I = 0.*pA
N_input.Idrive = 0.*mA
else:
run(dt_var*ms)
neuron.main.I = 1000.*pA
N_input.Idrive[nrIn+ttt] = 2000.*pA
run(3*ms-dt_var*ms)
neuron2.I = 0.*pA
N_input.Idrive[ttt] = 0.*mA
run(dt_var*ms)
neuron.I = 0.*pA
N_input.Idrive = 0.*mA
######
run(pair_interval*ms)
#store weight changes
weight_change1[jj] = 100.*(Syn_1.wampa[ttt] + 15.*(Syn_1.wampa[ttt]-init_weight))/init_weight
weight_change2[jj] = 100.*(Syn_2.wampa[ttt] + 15.*(Syn_2.wampa[ttt]-init_weight))/init_weight
run(5*ms)
print('Simulation finished!')
#####################################################
# Plots
#####################################################
#
titlestr = 'Data/'+morph[12:-8]+'_axonH_stdp_'+str_var+'_'+str(ttt)
# data1 = open(titlestr+'_AMPA_w1.txt','w')
# data2 = open(titlestr+'_AMPA_w2.txt','w')
# json.dump(weight_change1.tolist(),data1)
# json.dump(weight_change2.tolist(),data2)
# data1.close()
# data2.close()
x_values = np.array(list(reversed(-1.*dt_array)))
x_values = np.append(x_values,dt_array)
if paramS==0:
stitle = 'Prox Ca'
scolor = 'b'
clvar = [.5, .5, 1];
clvar2 = [.25, .25, 1];
else:
stitle = 'Dist Ca'
scolor = 'r'
clvar = [1, .5, .5];
clvar2 = [1, .25, .25];
fig = figure(figsize=(8, 5))
plt.plot(x_values,np.append(np.array(list(reversed(weight_change2))),weight_change1),'-',linewidth=2,color=scolor)
xlabel('Spike interval [ms]',fontsize=22)
ylabel('Normalised Weight [%]',fontsize=22)
plt.subplots_adjust(bottom=0.2,left=0.15,right=0.95,top=0.85)
title(stitle)
ylim([50,250])
# plt.savefig('./'+str(str_var)+'_june.eps', format='eps', dpi=1000)