// GENERATES AN IO CURVE OF RESPONSES TO BRIEF HIGH FREQUENCY STIMULATION 
// OF INCREASING NUMBER OF EXCITATORY INPUTS IN SLM
// TESTS FOR EFFECT OF GABAA ON INDUCTION OF NMDAR-MEDIATED DENDRITIC SPIKES 
//
// This script:
// - runs a family of two simulations: 
//		stimulation in stratum rediatum and stratum lacunosum moleculare,
//		same stimulation in the presence of AP5.
// For each simulation, voltages are displayed for the soma, a proximal tuft
// dendrite, and a distal tuft dendrite. Additionally, shape plots are used
// to illustrate the location of excitatory and inhibitory synapses.

// Numerical parameters for simulation.
tstop = 120 //  stop time of simulation
stimStart=20
dtime=0.2
stim_no=3
stim_inter=5
stim_delay=0
stim_noise=0
simul_iter=10

scale_factor1=0.2
scale_factor0=0.4

//seed of random number generator
randShift_inh = 2.357 
randShift_exc = 0.357 

//parameters for rectifying and non rectifying part of GABAsynapses
GABAweight_total=0.001 
GABAweight0=(GABAweight_total/5) 
GABAweight1=(4*GABAweight_total/5) 
GABAtau2=30
GABAtau1=1.0
GABArelP=1 
slope_factor=3
V50=-52

//establish rectifying and non-rectifying parts of GABA synapses
for ii=1,totVgatAt  {
    synGABA[ii-1].tau1=  GABAtau1/2
    synGABA[ii-1].tau2 =  GABAtau2/2
    synGABA[ii-1].slope_factor=0.3
    synGABA[ii-1].V50=-150

    synGABArect[ii-1].tau1=  GABAtau1
    synGABArect[ii-1].tau2 = GABAtau2
    synGABArect[ii-1].slope_factor=slope_factor
    synGABArect[ii-1].V50=V50
}

load_file("Generate_Stimulator.hoc")
load_file("Connect_Stimulator2ExcSyn.hoc")
load_file("Connect_Stimulator2InhSyn.hoc")

//setup recording vectors etc.
objref recv_tuft1, recv_tuft2, recv_tuft3, recv_soma, rect
objref  recv_tuft4, recv_tuft5, recv_trunk, recv_trunk2
strdef ColLabel, file_name1, file_name2, file_name3, file_name4
objref outfile, storeM_control, storeM_fast, storeM_noRect, storeM_fast_NR
load_file("Print-to-File.hoc")

ColLabel = "Membrane potential"
file_name1="C:/../Control.atf"
file_name3="C:/../fast_scaled.atf"
file_name2="C:/Daten/Modeling/Bloss_Spruston_Neuron_2016/CA1_pyr_JMS_v5/output/Syn_IO_GABAkinetics/nextround2/no_rect.atf"

recv_soma = new Vector()
recv_tuft1 = new Vector()
recv_tuft2 = new Vector()
recv_tuft3 = new Vector()
recv_tuft4 = new Vector()
recv_tuft5 = new Vector()
recv_trunk = new Vector()
recv_trunk2 = new Vector()
rect = new Vector()
recv_soma.record(&soma.sec.v(0.5),dtime)
rect.record(&t)

rec_conditions=8


storeM_fast = new Matrix(tstop/dtime+1,rec_conditions*simul_iter)
storeM_fast_NR = new Matrix(tstop/dtime+1,rec_conditions*simul_iter)
storeM_control = new Matrix(tstop/dtime+1,rec_conditions*simul_iter)
storeM_noRect = new Matrix(tstop/dtime+1,rec_conditions*simul_iter)

nExcAct_SLM = 150// number of excitatory synapses to activate, ~10% of 1464 total with spine density of 0.5/um2
nInhAct_SLM =  10 // number of inhibitory synapses to activate; ~3.5% of 280 total
nExcAct_SR = 0//200//number of excitatory synapses to activate, ~3% of 6000 total (oblique plus trunk; instead of 4920)
nInhAct_SR = 0//30 //number of inhibitory synapses to activate; 21% of 145 total 
// i.e. same proportion of inhibitory and excitatory synapses are activated at step 7

load_file("update_Synapses.hoc")

// Parameters for visualisation. dendInd1 and dendInd2 refer to dendritic
// indices that will have voltage traces plotted.
dendInd1 = 112 // left apical, SLM, with inhibition
dendInd2 = 149 // left apical, terminal neurite, SLM, no inhibition
dendInd3 = 179 // SLM, no inhibition
dendInd4 = 180 // SLM, no inhibition
dendInd5 = 176 // SLM, 2 GABA synapses
dendInd6 = 161 // apical trunk, just at SR/SLM border
dendInd7 = 85 //apical trunk, just at SR/SLM border
objref dend1Ref,dend2Ref,dend3Ref,dend4Ref,dend5Ref,dend6Ref,dend7Ref
	Cell[0].dend[dendInd1] {dend1Ref = new SectionRef()} 
	Cell[0].dend[dendInd2] {dend2Ref = new SectionRef()} 
	Cell[0].dend[dendInd3] {dend3Ref = new SectionRef()} 
	Cell[0].dend[dendInd4] {dend4Ref = new SectionRef()} 
	Cell[0].dend[dendInd5] {dend5Ref = new SectionRef()} 
	Cell[0].dend[dendInd6] {dend6Ref = new SectionRef()} 
	Cell[0].dend[dendInd7] {dend7Ref = new SectionRef()} 
    
recv_tuft1.record(&dend1Ref.sec.v(0.5),dtime)    
recv_tuft2.record(&dend2Ref.sec.v(0.5),dtime)    
recv_tuft3.record(&dend3Ref.sec.v(0.5),dtime)    
recv_tuft4.record(&dend4Ref.sec.v(0.5),dtime)    
recv_tuft5.record(&dend5Ref.sec.v(0.5),dtime)    
recv_trunk.record(&dend7Ref.sec.v(0.5),dtime)    
recv_trunk2.record(&dend6Ref.sec.v(0.5),dtime)

// Determine number of segments to track.
totSegs = 0
forall { for (x,0) { totSegs +=1 } }

// Create graphs for visualisation.
	objref voltBL_d1,voltBL_d2, voltAMPA 
	objref shapeExc[simul_iter/2],shapeInh

	voltBL_d1 = new Graph()
	voltBL_d1.addvar("soma.sec.v(0.5)",2,1)
	voltBL_d1.addvar("dend7Ref.sec.v(0.5)",3,1) 
	voltBL_d1.addvar("dend2Ref.sec.v(0.5)",4,1) 
	voltBL_d1.label("Non-linear")
	voltBL_d1.exec_menu("Keep Lines")
	voltBL_d1.size(stimStart-20,tstop,-75,-5)
    
	voltBL_d2 = new Graph()
	voltBL_d2.addvar("soma.sec.v(0.5)",2,1)
	voltBL_d2.addvar("dend7Ref.sec.v(0.5)",3,1) 
	voltBL_d2.addvar("dend1Ref.sec.v(0.5)",4,1) 
	voltBL_d2.label("Non-linear")
	voltBL_d2.exec_menu("Keep Lines")
	voltBL_d2.size(stimStart-20,tstop,-75,-5)
    
    strdef shape_label
	for (ii=1;ii<simul_iter/2+1; ii=ii+1){ //
        sprint(shape_label,"Exc, %d",2*ii)
        print shape_label
        shapeExc[ii-1] = new Shape()
        shapeExc[ii-1].label(shape_label)
    }
	shapeInh = new Shape()
	shapeInh.label("Inh")
    
// Create vectors that will store indices of active excitatory and 
// inhibitory synapses.
objref curExc_SLM, curExc_SR, curInh_SLM, curInh_SR  

// Vectors of normalised release probability 
objref norm_Pr_exc, norm_Pr_inh_SLM, ActSyn, ActSyn_inh

// Initialise simulations by adding channels and assigning genotypes to
// inhibitory synapses.
initChannels()
seedGenotypes()
	
norm_Pr_exc = new Vector(5) //, [1, 1.55,	1.83,	1.93,	2.03])
ActSyn = new Vector(5)
norm_Pr_exc.x[0] = 1
norm_Pr_exc.x[1] = 1.5 
norm_Pr_exc.x[2] = 1.8
norm_Pr_exc.x[3] = 1.9 
norm_Pr_exc.x[4] = 2.0

// DO SIMULATION AND PLOT RESULTS.
// turn off cvode variable time step in order to avoid errors during repetitive simulations
cvode_active(0)


// run simulation with inhibition
curGr = graphList[0].append(voltBL_d1)
curGr = graphList[0].append(voltBL_d2)
norm_Pr_inh_SLM = new Vector(5) //
ActSyn_inh = new Vector(5)
norm_Pr_inh_SLM.x[0] = 1
norm_Pr_inh_SLM.x[1] = 1
norm_Pr_inh_SLM.x[2] = 1   
norm_Pr_inh_SLM.x[3] = 1 
norm_Pr_inh_SLM.x[4] = 1

activateInhibition(cellList,-1,1) // clear inhibition
curInh_SLM = activateInhibition(tuftList,nInhAct_SLM,randShift_inh,"vgat",1) 
ActSyn_inh = set_InhSyn_syn_fixed(curInh_SLM, randShift_inh/4, nInhAct_SLM, shapeInh) 
print "In SLM, the number of activated inhibitory synapses at each pulse are: ",ActSyn_inh.x[0], ", " , ActSyn_inh.x[1], ", " , ActSyn_inh.x[2], ", " , ActSyn_inh.x[3], ", " , ActSyn_inh.x[4]

for(jj=1; jj<simul_iter+1; jj=jj+1){
    print jj
	// Clear excitation and then turn on select synapses.
	activateExcitation(cellList,-1,1) // clear excitation
	curExc_SLM = activateExcitation(tuftList,jj*nExcAct_SLM,randShift_exc) // activate excitatory synapses
	
    shape_no=(jj/2)-1
    if (jj%2==1){shape_no=(jj-1)/2}
	
    ActSyn = set_gluSyn_fixed_N(curExc_SLM, randShift_exc/2, 0.1, norm_Pr_exc, shapeExc[shape_no])
	print "In SLM, the number of activated synapses at each pulse are: ",ActSyn.x[0], ", " , ActSyn.x[1], ", " , ActSyn.x[2], ", " , ActSyn.x[3], ", " , ActSyn.x[4]
	    
    init() 
    run()
    storeM_control.setcol(jj-1,recv_soma)
    storeM_control.setcol(jj-1+simul_iter,recv_tuft1) 
    storeM_control.setcol(jj-1+2*simul_iter,recv_tuft2) 
    storeM_control.setcol(jj-1+3*simul_iter,recv_tuft3)
    storeM_control.setcol(jj-1+4*simul_iter,recv_tuft4) 
    storeM_control.setcol(jj-1+5*simul_iter,recv_tuft5) 
    storeM_control.setcol(jj-1+6*simul_iter,recv_trunk) 
    storeM_control.setcol(jj-1+7*simul_iter,recv_trunk2)
        
    reset_gluSyn()
}
graphList[0].remove_all()


// SAVE OUTPUT
//print2file(storeM_control,file_name1,ColLabel)


// run simulation without rectification inhibition
load_file("tuft_NMDA_spike_noRect.hoc")

// run simulation with PV-like dynamics
load_file("tuft_NMDA_spike_fast.hoc")