// GENERATES AN IO CURVE OF RESPONSES TO STIMULATION OF INCREASING NUMBER OF EXCITATORY INPUTS IN SR AND SLM
// TESTS FOR CONTRIBUTION OF NMDARs AND EFFECT OF NONLINEAR GABAA
//
// 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 = 220 //  stop time of simulation
stimStart=20
dtime=0.2
stim_no=5
stim_inter=20
stim_delay=0
stim_noise=0
simul_iter=8//10

//seed of random number generator
randShift = 6.37


//parameters for rectifying and non rectifying part of GABAsynapses
GABAweight_total=0.0007
GABAweight0=GABAweight_total/5 
GABAweight1=4*GABAweight_total/5  
GABAtau2=30
GABAtau1=1.0
GABArelP=0.5
slope_factor=4 
V50=-60

//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_obl1, recv_obl2, recv_obl3
strdef ColLabel, file_name1, file_name2, file_name3
objref outfile, storeM_noInh, storeM_control, storeM_redRect
load_file("Print-to-File.hoc")

ColLabel = "Membrane potential"
file_name1="C:/../No_Inh.atf"
file_name2="C:/../Control.atf"
file_name3="C:/../reduced_rect.atf"

recv_soma = new Vector()
recv_tuft1 = new Vector()
recv_tuft2 = new Vector()
recv_tuft3 = new Vector()
recv_obl1 = new Vector()
recv_obl2 = new Vector()
recv_obl3 = new Vector()
rect = new Vector()
recv_soma.record(&soma.sec.v(0.5),dtime)
rect.record(&t)

rec_conditions=14

storeM_noInh = new Matrix(tstop/dtime+1,rec_conditions*simul_iter) 
storeM_control = new Matrix(tstop/dtime+1,rec_conditions*simul_iter)
storeM_redRect = new Matrix(tstop/dtime+1,rec_conditions*simul_iter)

nExcAct_SLM = 50// number of excitatory synapses to activate, ~3.3% of 1464 total with spine density of 0.5/um2
nInhAct_SLM =  95 // number of inhibitory synapses to activate; 34% of 280 total
nExcAct_SR = 200//number of excitatory synapses to activate, ~3% of 6000 total (oblique plus trunk; instead of 4920)
nInhAct_SR = 50 //number of inhibitory synapses to activate; 34% 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 = 107 // left apical, SLM
dendInd2 = 149 // left apical, terminal neurite, SLM
dendInd3 = 176 // SLM
dendInd4 = 186 // SR
dendInd5 = 144 // SR
dendInd6 = 136 // SR
objref dend1Ref,dend2Ref,dend3Ref,dend4Ref,dend5Ref,dend6Ref
	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()} 
    
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_obl1.record(&dend4Ref.sec.v(0.5),dtime)    
recv_obl2.record(&dend5Ref.sec.v(0.5),dtime)    
recv_obl3.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,voltAMPA 
	objref shapeExc[simul_iter],shapeInh

	voltBL = new Graph()
	voltBL.addvar("soma.sec.v(0.5)",2,1)
	voltBL.addvar("dend2Ref.sec.v(0.5)",4,1) 
	voltBL.label("No inhib")
	voltBL.exec_menu("Keep Lines")
	voltBL.size(stimStart-20,tstop,-75,-15)
	
	voltAMPA = new Graph()
	voltAMPA.addvar("soma.sec.v(0.5)",2,1)
	voltAMPA.addvar("dend2Ref.sec.v(0.5)",4,1)  
	voltAMPA.label("AMPA only")
	voltAMPA.exec_menu("Keep Lines")
	voltAMPA.size(stimStart-20,tstop,-75,-15)
    
    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 release of normalised release probability 
objref norm_Pr_exc, norm_Pr_inh_SR, norm_Pr_inh_SLM, ActSyn, ActSyn_inh, SynAct

// Initialise simulations by adding channels and assigning genotypes to
// inhibitory synapses.
initChannels()
seedGenotypes()

// Implement short-term plasticity at excitatory synapses
norm_Pr_exc = new Vector(5) //
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 without inhibition for different numbers of glu inputs
load_file("SynStim_SR_SLM_noInh.hoc")

// activate inhibitory inputs
activateInhibition(cellList,-1,1) // clear inhibition
curInh_SR = activateInhibition(radiatumList,nInhAct_SR,randShift,"vgat",1)
curInh_SLM = activateInhibition(tuftList,nInhAct_SLM,randShift,"vgat",1) 
   
// Implement short-term plasticity at inhibitory synapses
norm_Pr_inh_SR = new Vector(5) //
ActSyn_inh = new Vector(5)
norm_Pr_inh_SR.x[0] = 1
norm_Pr_inh_SR.x[1] = 0.7   
norm_Pr_inh_SR.x[2] = 0.7  
norm_Pr_inh_SR.x[3] = 0.7
norm_Pr_inh_SR.x[4] = 0.7

ActSyn_inh = set_InhSyn_fixed_N(curInh_SR, randShift/3, GABArelP, norm_Pr_inh_SR, shapeInh) 
print "In SR, 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]

// Implement short-term plasticity at inhibitory synapses
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] = 0.6 
norm_Pr_inh_SLM.x[2] = 0.5   
norm_Pr_inh_SLM.x[3] = 0.4 
norm_Pr_inh_SLM.x[4] = 0.4

ActSyn_inh = set_InhSyn_fixed_N(curInh_SLM, randShift/4, GABArelP, norm_Pr_inh_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]

// run simulation with normal inhibition
load_file("SynStim_SR_SLM_control.hoc")
// run simulation with reduced rectifying inhibition
load_file("SynStim_SR_SLM_redInh.hoc")