/*----------------------------------------------------------------------------
	Andrew Knox 2014
	Newer synapses

----------------------------------------------------------------------------*/

nsynapses = 100
nstimuli = 10
nhfburst = 10 //10
timestimuli = 333
timehfburst = 3 //2.5

//----------------------------------------------------------------------------
//  load and define general graphical procedures
//----------------------------------------------------------------------------

load_file("nrngui.hoc")

objectvar g[20]			// max 20 graphs
objectvar g1, g2, g3, g4
ngraph = 0
objref PYtimevec, PYidvec, INtimevec, INidvec, REtimevec, REidvec, TCtimevec, TCidvec, recncs, recveclist

field = 0
gababi = 0
ampai = 0

proc addgraph() { local ii	// define subroutine to add a new graph
				// addgraph("variable", minvalue, maxvalue)
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph(0)
	g[ii].view(0,1,0,1, int(ii/2)*550+80, ii%2*450+80, 400, 300)
	g[ii].size(tstart,tstop,$2,$3)
	g[ii].xaxis()
	g[ii].yaxis()
	g[ii].addvar($s1,1,0)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
}

proc addtext() { local ii	// define subroutine to add a text graph
				// addtext("text")
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph()
	g[ii].size(0,tstop,0,1)
	g[ii].xaxis(3)
	g[ii].yaxis(3)
	g[ii].label(0.1,0.8,$s1)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
	text_id = ii
}

proc addline() {		// to add a comment to the text window
				// addline("text")
	g[text_id].label($s1)
}

if(ismenu==0) {
  nrnmainmenu()			// create main menu
  nrncontrolmenu()		// create control menu
  ismenu=1
}



objref membranedatafile
membranedatafile = new File()

//----------------------------------------------------------------------------
//  transient time
//----------------------------------------------------------------------------

trans = 00

print " "
print ">> Transient time of ",trans," ms"
print " "


//----------------------------------------------------------------------------
//  setup simulation parameters
//----------------------------------------------------------------------------

Dt = 0.1
npoints = 40000

objectvar Sim			// create vector of simulation points
Sim = new Vector(npoints)

dt = 0.1			// must be submultiple of Dt
tstart = trans
tstop = trans + npoints * Dt
runStopAt = tstop
steps_per_ms = 1/Dt
celsius = 36
v_init = -70.5			//originally 70




//----------------------------------------------------------------------------
//  create PY cell
//----------------------------------------------------------------------------

load_file("sPY.tem")		// read geometry file

objectvar PY		// create PY cells
objectvar PYVtrace


PY = new sPY()
PYVtrace = new Vector()
PYVtrace.record(&PY.soma.v(0.5))


//----------------------------------------------------------------------------
//  setup stimulus vectors
//----------------------------------------------------------------------------

//load_file("vecevent.mod")



objectvar gabaAtimevec[nsynapses]
objectvar gabaBtimevec[nsynapses]
objectvar ampatimevec[nsynapses]
objectvar nmdatimevec[nsynapses]

objectvar nmdavec[nsynapses]

objref VgabaA[nsynapses] 
objref VgabaB[nsynapses] 
objref Vampa[nsynapses] 

double Vnmda[nsynapses] 

objectvar pygabaB[nsynapses]
objectvar pynmda[nsynapses]

objectvar jitter
jitter = new Random()
jitter.uniform(-1,1)

for ctr=0,nsynapses-1 {
   gabaAtimevec[ctr] = new Vector()
   gabaBtimevec[ctr] = new Vector()
   ampatimevec[ctr] = new Vector()
   nmdatimevec[ctr] = new Vector()

   nmdavec[ctr] = new Vector()

   VgabaA[ctr] = new VecStim()
   VgabaB[ctr] = new VecStim()
   Vampa[ctr] = new VecStim()
   Vnmda[ctr] = 0

   nmdavec[ctr].append(-60)
   nmdatimevec[ctr].append(0)

   for i=1, nstimuli-1 {
     for j=1, nhfburst-1 {
       stime = i*timestimuli+j*timehfburst+jitter.repick()
       //stime = i*timestimuli+j*timehfburst
       gabaAtimevec[ctr].append(stime)
       
       stime = i*timestimuli+j*timehfburst+jitter.repick()
       //stime = i*timestimuli+j*timehfburst
       gabaBtimevec[ctr].append(stime)

       stime = i*timestimuli+j*timehfburst+jitter.repick()
       //stime = i*timestimuli+j*timehfburst
       ampatimevec[ctr].append(stime)

       stime = i*timestimuli+j*timehfburst+jitter.repick()
       //stime = i*timestimuli+j*timehfburst
       nmdatimevec[ctr].append(stime)
       nmdavec[ctr].append(10)
       nmdatimevec[ctr].append(stime+2*Dt)
       nmdavec[ctr].append(-60)       
     }
   }
   nmdavec[ctr].play(&Vnmda[ctr],nmdatimevec[ctr])
}



//----------------------------------------------------------------------------
//  setup synapses 
//----------------------------------------------------------------------------

objectvar pygabaA
objectvar pyampa

objref netcon

//   pygabaA = new GABAa_S()
//   pyampa = new AMPA_S()
//   pynmda = new NMDA() 


//   PY.soma pygabaA.loc(0.5)
//   PY.soma pyampa.loc(0.5)
//   PY.soma pynmda.loc(0.5)

   for i=0,nsynapses-1 {
      pygabaB[i] = new GABAb_S()
      pynmda[i] = new NMDA()
      PY.soma pygabaB[i].loc(0.5)
      PY.soma pynmda[i].loc(0.5)
      PY.gababpost.append(pygabaB[i])
   }	


Alpha_GABAa_S = 20 		// from diffusion model
Beta_GABAa_S = 0.16 
Cmax_GABAa_S = 0.5 		// short pulses
Cdur_GABAa_S = 0.3 
Erev_GABAa_S = -80 		// Rinzel's Erev


K1_GABAb_S	= 0.09	//0.52		(/ms mM) forward binding rate to receptor
K2_GABAb_S	= 0.0012//0.0045	(/ms)	backward (unbinding) rate of receptor
K3_GABAb_S	= 0.18 	//		(/ms)	rate of G-protein production
K4_GABAb_S	= 0.034	//		(/ms)	rate of G-protein decay
KD_GABAb_S	= 100	//		dissociation constant of K+ channel
n_GABAb_S	= 4	//		nb of binding sites of G-protein on K+
Erev_GABAb_S	= -95 	//	(mV)	reversal potential (E_K)

Cmax_GABAb_S = 0.5 	// short pulses
Cdur_GABAb_S = 0.3 


Alpha_AMPA_S = 0.94		// kinetics from simplex with short pulses
Beta_AMPA_S = 0.18
Cmax_AMPA_S = 0.5
Cdur_AMPA_S = 0.3
Erev_AMPA_S = 0

Cmax_NMDA = 0.5	//1	//(mM)		: max transmitter concentration
Cdur_NMDA = 0.3 //1		//(ms)		: transmitter duration (rising phase)
Alpha_NMDA = 0.11 //.072	//(/ms mM)	: forward (binding) rate
Beta_NMDA = 0.0066 	//(/ms)		: backward (unbinding) rate
Erev_NMDA = 0		//(mV)
Mg_NMDA = 2		//(mM)		: external magnesium

proc assign_synapses() {	// procedure to assign syn conductances 
				// params: 1=gabaA, 2=gabaB, 3=ampa, 4=nmda
   for i=0,nsynapses-1 {  
	PY.INgabaalist.append(new NetCon(VgabaA[i], PY.gabaapost, 0, 0, $1))
	PY.INgabablist.append(new NetCon(VgabaB[i], PY.gababpost.object(i), 0, 0, 1))
	PY.gababpost.object(i).gmax = $2
	PY.PYlist.append(new NetCon(Vampa[i], PY.ampapostTC, 0, 0, $3))
        setpointer pynmda[i].pre, Vnmda[i]
	pynmda[i].gmax = $4   
   }
}

//assign_synapses(0.0015,0.004,0.004)
assign_synapses(0.0015,0.004,0.004,0.001)
//assign_synapses(0.00,0.004,0.005,0.00)
//assign_synapses(1.5,4,4,1)
//assign_synapses(0,0.004,0.005)

for i=0,nsynapses-1 { 
   VgabaA[i].play(gabaAtimevec[i])
   VgabaB[i].play(gabaBtimevec[i])
   Vampa[i].play(ampatimevec[i])
}



//----------------------------------------------------------------------------
//  add graphs
//----------------------------------------------------------------------------

//addgraph("tcB[0][0].g",0,0.05)

//addgraph("TC[0].soma.o2_iar",0,1)
//addgraph("TC[0].soma.p1_iar",0,1)
//addgraph("TC[15].it.h",0,1)
//addgraph("TC[15].it.h_inf",0,1)


strdef gtxt


sprint(gtxt,"PY.soma.v(0.5)")
addgraph(gtxt,-90,10)



addgraph("field",-5,2)

//addgraph("VgabaB[0]",-70,15)
//addgraph("Vampa[1]",-70,15)
//addgraph("Vampa[2]",-70,15)

addgraph("ampai",-20,0.1)
//addgraph("PY.ampapost.g",-0.1,0.1)
//addgraph("PY.ampapost.Ron",-0.1,0.1)
//addgraph("PY.ampapost.Roff",-0.1,0.1)
//addgraph("PY.gabaapost.i",-0.1,0.1)
//addgraph("PY.gababpost.object(0).i",0,0.005)
addgraph("gababi",0,5)
//addgraph("PY.gababpost.R",-0.1,0.4)
//addgraph("PY.gababpost.Gn",-0.1,0.4)
//addgraph("PY.gababpost.Roff",-0.1,0.4)
//addgraph("PY.INgabablist.object(0).weight[1]",-0.001,0.001)
//addgraph("PY.INgabablist.object(0).weight[2]",325,375)
//addgraph("PY.ampapost.synon",-0.001,0.001)


//-----------------------------------------------------------------------------
// Write recorded membrane data to file
//-----------------------------------------------------------------------------

proc writedatafile()  {

  membranedatafile.wopen("membrane_data.txt")
  PYVtrace.printf(membranedatafile) 
  membranedatafile.close()
}

////////////////////////////////////////////////////////////////////////
//  Code for dealing with field potentials
///////////////////////////////////////////////////////////////////////

func xfield()  {local i, j, tmp, Re, x, cw   
//$1 is Re, $2 is distance from nearest neuron
   total_field = 0
   tmp = 0
   Re = $1
   x = $2

   for i=0, nsynapses-1  {
	tmp += pygabaB[i].i
	tmp += pynmda[i].g
   }

   return (PY.ampapostTC.i + PY.gabaapost.i + tmp) * Re/4/PI/x
//   return (PY.ampapostTC.i + PY.gabaapost.i + PY.soma.ik_im + tmp) * Re/4/PI/x
}

func gababcur()  {local i, tmp
   tmp = 0
   Re = $1
   x = $2
   for i=0, nsynapses-1 {
      tmp += PY.gababpost.object(i).i
   }
   return tmp * Re/4/PI/x
}

func ampacur()  {
   Re = $1
   x = $2
   
   return PY.ampapostTC.i * Re/4/PI/x
}

proc advance() {
   fadvance()
   field = xfield(230,5)
   gababi = gababcur(230,5)
   ampai = ampacur(230,5)
}

proc printweight() {
   print PY.INgabablist.object(0).weight
   print PY.INgabablist.object(0).weight[0]
   print PY.INgabablist.object(0).weight[1]
   print PY.gababpost.K2
}