// GENERATE A FAMILY OF SINGLE AND BURST PSPS DEPENDING ON E-I RATIO
// JM Schulz, 2020


objref cvode
cvode = new CVode()
cvode.active(1)


strdef ColLabel, file_name1, file_name2
objref shapeSyn 
//save output in python vectors
nrnpython("import numpy as N")
nrnpython("import matplotlib.pyplot as plt")
nrnpython("from mpl_toolkits import mplot3d")

objref p
p = new PythonObject()
nrnpython("Peak_PSP=N.zeros((0,))")
nrnpython("Peak_PSP_AP5=N.zeros((0,))")
nrnpython("X=N.zeros((0,))")
nrnpython("Y=N.zeros((0,))")
nrnpython("Arr_Peak_PSP=N.zeros((21,61))")
nrnpython("Arr_Peak_PSP_AP5=N.zeros((21,61))")
nrnpython("Arr_X=N.zeros((21,61))")
nrnpython("Arr_Y=N.zeros((21,61))")
nrnpython("x=0")
nrnpython("y=0")
nrnpython("peak=0")

load_file("Calc_Func.hoc")

// Numerical parameters for simulation.
tstop = 350 // stop time of simulation
stimStart=20
dtime=0.2
gluSyn_iter=21
GabaSyn_iter=61

// GABA reversal
inhRev = -35

ampaWeight=0.00002 // 0.00002 = 1/7 of mature synapse in CA1 (Schulz et al. 2018); Li et al., 2017: 2 pA = 0.025 nS
nmdaWeight=7*ampaWeight // Li et al., 2017
GABAweight_total=0.00024

alpha_vspom=-0.062
v0_block=10

//synGABA[0] is not rectifying, synGABA[1] is rectifying
GABAweight0= GABAweight_total/5 
GABAweight1= 4*GABAweight_total/5  
GABAtau2=55
GABAtau1=3

shapeSyn = new Shape()
shapeSyn.label("Syn")

for (ii=0; ii<n_GABA_syn; ii=ii+2){
    synGABA[ii].tau1= 0.5
    synGABA[ii].tau2 = 10 
    ncGABA[ii].weight = GABAweight0 
    synGABA[ii].slope_factor=0.3
    synGABA[ii].V50=-150
    shapeSyn.point_mark(synGABA[ii],3,"O",6)

    synGABA[ii+1].tau1= GABAtau1
    synGABA[ii+1].tau2 = GABAtau2
    ncGABA[ii+1].weight = GABAweight1
	if (GABA_option==1) {synGABA[ii+1].slope_factor=7    //Analysis of VC data in Schulz et al., 2019
    synGABA[ii+1].V50=  -45}
	if (GABA_option==2) {synGABA[ii+1].slope_factor=0.3 
    synGABA[ii+1].V50=  250}
}

for (ii=0; ii<n_glu_syn; ii=ii+1){
    ncAMPA[ii].weight = ampaWeight
    ncNMDA[ii].weight = nmdaWeight
    nsAMPA[ii].number = 3
    nsAMPA[ii].start = stimStart
    nsNMDA[ii].number = 3
    nsNMDA[ii].start = stimStart
    shapeSyn.point_mark(synAMPA[ii],2,"O",6)
}

//setup recording vectors
objref recv, rect, reci_NMDA
recv = new Vector()
rect = new Vector()

access soma[1]
recv.record(&soma[1].v(0.5),dtime)
rect.record(&t)

// Create graphs for visualisation.
objref Inh_intact_5, Inh_intact_10, Inh_intact_15, Inh_intact_20, NMDA_plus_AMPA
objref Inh_only

Inh_intact_5 = new Graph() 
Inh_intact_5.label("5 GABA syn")
Inh_intact_5.addvar("soma[1].v(0.5)",2,1)
Inh_intact_5.exec_menu("Keep Lines")
Inh_intact_5.size(0,tstop,-85,-15)

Inh_intact_10= new Graph() 
Inh_intact_10.label("10 GABA syn")
Inh_intact_10.addvar("soma[1].v(0.5)",2,1)
Inh_intact_10.exec_menu("Keep Lines")
Inh_intact_10.size(0,tstop,-85,-15)

Inh_intact_15 = new Graph() 
Inh_intact_15.label("30 GABA syn")
Inh_intact_15.addvar("soma[1].v(0.5)",2,1)
Inh_intact_15.exec_menu("Keep Lines")
Inh_intact_15.size(0,tstop,-85,-15)

Inh_intact_20 = new Graph() 
Inh_intact_20.label("20 GABA syn")
Inh_intact_20.addvar("soma[1].v(0.5)",2,1)
Inh_intact_20.exec_menu("Keep Lines")
Inh_intact_20.size(0,tstop,-85,-15)

NMDA_plus_AMPA = new Graph()
NMDA_plus_AMPA.label("gabazine")
NMDA_plus_AMPA.addvar("soma[1].v(0.5)",2,1)
NMDA_plus_AMPA.exec_menu("Keep Lines")
NMDA_plus_AMPA.size(0,tstop,-85,-15)

Inh_only = new Graph()
Inh_only.label("AP5 + NBQX")
Inh_only.addvar("soma[1].v(0.5)",2,1)
Inh_only.exec_menu("Keep Lines")
Inh_only.size(0,tstop,-85,-15)

// Adjust GABA synapse dynamics
for (ii=0; ii<n_GABA_syn; ii=ii+1){
    synGABA[ii].isOn=0
    synGABA[ii].e = inhRev

    nsGABA[ii].interval = 20
    nsGABA[ii].number = 3
    nsGABA[ii].start = stimStart
    nsGABA[ii].noise = 0
    ncGABA[ii].delay = 0
}

for (ii=0; ii<n_glu_syn; ii=ii+1){
    synAMPA[ii].isOn=0
    synNMDA[ii].isOn=0
    synAMPA[ii].tau1= 0.2  // Jarsky et al., 2005
    synAMPA[ii].tau2 = 2 // Jarsky et al., 2005

    nsAMPA[ii].interval = 20
    nsAMPA[ii].noise = 0
    ncAMPA[ii].delay = 0

    synNMDA[ii].tau1= 3  
    synNMDA[ii].tau2 = 35
    
    nsNMDA[ii].interval = 20
    nsNMDA[ii].noise = 0
    ncNMDA[ii].delay = 0
    synNMDA[ii].alpha_vspom=alpha_vspom
    synNMDA[ii].v0_block=v0_block
}

simul_no=0    
// run simulation for increasing synapse numbers
for(gg=0; gg<GabaSyn_iter; gg=gg+1){
    if(gg>0){
        synGABA[2*(gg-1)].isOn=1
        synGABA[2*(gg-1)+1].isOn=1
    }
    for(jj=0; jj<gluSyn_iter; jj=jj+1){
        synAMPA[jj-1].isOn=0
        synNMDA[jj-1].isOn=0
    }
    for(jj=0; jj<gluSyn_iter; jj=jj+1){
        if(jj>0){
            synAMPA[jj-1].isOn=1
            synNMDA[jj-1].isOn=1
        }      
        
        if (gg==5 && jj%4==1){curGr = graphList[0].append(Inh_intact_5)}
        if (gg==10 && jj%4==1){curGr = graphList[0].append(Inh_intact_10)}
        if (gg==30 && jj%4==1){curGr = graphList[0].append(Inh_intact_15)}
        if (gg==20 && jj%4==1){curGr = graphList[0].append(Inh_intact_20)}
        if (gg==0 && jj%4==1){curGr = graphList[0].append(NMDA_plus_AMPA) }
        if (jj==0 && gg%5==1){curGr = graphList[0].append(Inh_only) }
		
        finitialize()
        run()
		
        if (gg==5 && jj%4==1){graphList[0].remove(curGr-1)}
        if (gg==10 && jj%4==1){graphList[0].remove(curGr-1)}
        if (gg==30 && jj%4==1){graphList[0].remove(curGr-1)}
        if (gg==20 && jj%4==1){graphList[0].remove(curGr-1)}
        if (gg==0 && jj%4==1){graphList[0].remove(curGr-1)}
        if (jj==0 && gg%5==1){graphList[0].remove(curGr-1)}
		
        p.x=jj
        p.y=gg
        p.peak=get_max(recv,int(stimStart/dtime),int(tstop/dtime)-1)+80
        nrnpython("Peak_PSP=N.append(Peak_PSP,peak)")
        nrnpython("Arr_Peak_PSP[int(x),int(y)]=peak")
        
        nrnpython("X=N.append(X,x)")
        nrnpython("Y=N.append(Y,y)")
        
        nrnpython("Arr_X[int(x),int(y)]=x")
        nrnpython("Arr_Y[int(x),int(y)]=y")
		
    }       
}
print "simulation is done"

nrnpython("fig1 = plt.figure()")
nrnpython("ax = plt.axes(projection='3d')")
nrnpython("ax.plot_surface(Arr_X, Arr_Y, Arr_Peak_PSP, rstride=1, cstride=1, cmap='jet', edgecolor='none')")
nrnpython("ax.scatter(X,Y, Peak_PSP)")
nrnpython("ax.set_xlim([21,-1])")
nrnpython("ax.set_ylim([-0.5,61])")
nrnpython("ax.set_xlabel('# Glu synapses')")
nrnpython("ax.set_ylabel('# GABA synapses')")
nrnpython("ax.set_zlabel('Peak Depol (mV)')")

nrnpython("plt.show()")