// Firing analysis.hoc analyzes model firing generated by firing simulator.hoc.
// Number of APs, raster plot, interspike interval(ISI) and 1stAP latency against
// current intensity and currents normalized at threshold current (1T) are measured.
// Threshold current, 1stAP latency, number of APs at 1.5T and 2T, ISI(last)/ISI(1st) ratio
// are printed in command prompt area, which are used for comprehensive analysis (Figures 7A-E). 
//
// Takaki Watanabe
// wtakaki@m.u-tokyo.ac.jp

objref vdat, v1dat, v2dat, v3dat, tdat, sdat
objref spikes,spiket,spikes2,thr1, thr1h, thr2, thr2h, thr3, thr2000, thr3000, tdatfile
objref thresholdcur
objref gnadat, gkcnadat, gkcnab2dat, gkhtdat, gkadat, gkcnqdat
objref datgna, datgkcna, datgkcnab2, datgkht, datgka, datgkcnq

psize=10000
vdat = new Vector(psize,0)
v1dat = new Vector(psize,0)
v2dat = new Vector(psize,0)
v3dat = new Vector(psize,0)
tdat = new Vector(psize,0)
sdat = new Vector(psize,0)
spikes = new Vector(psize,0)
spikes2 = new Vector(psize,0)
spiket = new Vector(psize,0)
gnadat = new Vector(psize,0)
gkcnadat = new Vector(psize,0) 
gkcnab2dat = new Vector(psize,0) 
gkhtdat  = new Vector(psize,0)
gkadat = new Vector(psize,0) 
gkcnqdat = new Vector(psize,0)

set_original()
vdat.record(&soma.v(0.5))
v1dat.record(&soma.i_cap(0.5))
tdat.record(&t)
gnadat.record(&soma.gna_na(0.5))
gkcnadat.record(&soma.gkcna_kcna(0.5))
gkcnab2dat.record(&soma.gkcnab2_kcnab2(0.5))
gkhtdat.record(&soma.gkht_kht(0.5))
gkadat.record(&soma.gka_ka(0.5))
gkcnqdat.record(&soma.gkcnq_kcnq(0.5))

thr1 = new File()
thr1h = new File()
thr2 = new File()
thr2h = new File()
thr3 = new File()
thr2000 = new File()
thr3000 = new File()
tdatfile = new File()
thresholdcur = new File()
datgna = new File()
datgkcna = new File()
datgkcnab2 = new File()
datgkht = new File() 
datgka = new File() 
datgkcnq = new File()

yt.x[4]=0.2
yb.x[4]=-0.2
yt.x[5]=2
yb.x[5]=-0.2


proc run3() {
stoprun=0
	ana[4].beginline(color5,1)
	ana[5].beginline(color5,1)
	ana[8].beginline(color5,1)
	ana[9].beginline(color5,1)
v2dat=v1dat.c.deriv()
v3dat=v2dat.c.deriv()

sdat = vdat.c
spikes.spikebin(sdat,-20)
spiket.indvwhere(spikes,">", 0)
spikes2 = spikes.c.indvwhere(spikes,">",0)
spiken = spiket.size 

ana[3].mark(ic.amp,spiken,"O",3,color5,1)
ana[3].line(ic.amp,spiken)


if(spiken!=0&&threcur==0.01){
threcur=icur.x[i]
print threcur // threshold current
wopen("thresholdcur.dat")
fprint("%d\n",threcur)
wopen()
normcur=1
}

if(gkcnqbar_kcnq1==0){
threcur0=threcur
}

if(icur.x[i]>=threcur*1.48 && icur.x[i]<= threcur*1.52 && normcur==1 && drwn2==0){
normcur=1.5
}else if(icur.x[i]==threcur*2 && drwn3==0){
normcur=2
}else if(icur.x[i]>=threcur*2.48 && icur.x[i]<= threcur*2.52 && normcur==2 && drwn4==0){
normcur=2.5
}else if(icur.x[i]==threcur*3 && drwn5==0){
normcur=3
}

if(spiken>=1){
for k=0,spiken-1 {
   ana[4].mark(spikes2.c.x[k]*dt,ic.amp,"O",3,color5,1)
   ana[4].line(spikes2.c.x[k]*dt,ic.amp)
   }
	if(spiken>1){
for k=1,spiken-1 {
	ana[5].mark(k+1,spikes2.x[k]*dt-spikes2.x[k-1]*dt,"O",3,color5,1)
	ana[5].line(k+1,spikes2.x[k]*dt-spikes2.x[k-1]*dt)
	}
	}
	ana[6].mark(ic.amp,spikes2.x[0]*dt-icdel,"O",3,color5,1)
	ana[6].line(ic.amp,spikes2.x[0]*dt-icdel)
}

if (normalize==1 && normcur>=1){
ana[7].mark(normcur,spiken,"O",3,color5,1)
ana[7].line(normcur,spiken)
	if(normcur==1.5&&norm2==0||normcur==2&&norm2==1){
	print spiken // number of APs at 1.5T and 2T
norm2=norm2+1
}

for k=0,spiken-1 {
   ana[8].mark(spikes2.c.x[k]*dt,normcur,"O",3,color5,1)
   ana[8].line(spikes2.c.x[k]*dt,normcur)
   }
	if(spiken>1){
for k=1,spiken-1 {
	ana[9].mark(k+1,spikes2.x[k]*dt-spikes2.x[k-1]*dt,"O",3,color5,1)
	ana[9].line(k+1,spikes2.x[k]*dt-spikes2.x[k-1]*dt)
	}
	}
	ana[10].mark(normcur,spikes2.x[0]*dt-icdel,"O",3,color5,1)
	ana[10].line(normcur,spikes2.x[0]*dt-icdel)
}else if (normalize==0 && icur.x[i]==threcur && drwn1==0||(icur.x[i]>=threcur*1.48 && icur.x[i]<= threcur*1.52 && drwn2==0)||icur.x[i]==threcur*2 && drwn3==0||(icur.x[i]>=threcur*2.48 && icur.x[i]<= threcur*2.52 && drwn4==0)||icur.x[i]==threcur*3 && drwn5==0){

ana[7].mark(normcur,spiken,"O",3,color5,1)
ana[7].line(normcur,spiken)

for k=0,spiken-1 {
   ana[8].mark(spikes2.c.x[k]*dt,normcur,"O",3,color5,1)
   ana[8].line(spikes2.c.x[k]*dt,normcur)
   }
	if(spiken>1){
for k=1,spiken-1 {
	ana[9].mark(k+1,spikes2.x[k]*dt-spikes2.x[k-1]*dt,"O",3,color5,1)
	ana[9].line(k+1,spikes2.x[k]*dt-spikes2.x[k-1]*dt)
	}
	}
	ana[10].mark(normcur,spikes2.x[0]*dt-icdel,"O",3,color5,1)
	ana[10].line(normcur,spikes2.x[0]*dt-icdel)
}

    ana[3].size(iint, ic.amp*1.2, yb.x[6], yt.x[6])
	ana[4].size(-xscale1*0.1, xscale1,-ic.amp*0.2, ic.amp*1.2)
	ana[5].size(0, spiken, yb.x[8], yt.x[8])
	ana[6].size(-ic.amp*0.1, ic.amp*1.2, yb.x[9], yt.x[9])
	ana[7].size(1, 3, yb.x[6], yt.x[6])
	ana[8].size(-xscale1*0.1, xscale1,0, 4)
	ana[9].size(0, spiken, yb.x[8], yt.x[8])
	ana[10].size(1, 3, yb.x[9], yt.x[9])
	
	for p=3,10{
	ana[p].flush()
	   ana[p].exec_menu("View = plot")
	ana[p].exec_menu("10% Zoom out")
	}


if(normcur==1 && drwn1==0){
    vdat.c.line(ana[num],tdat,color5,1)
	ana[num].flush
	thr1.wopen("thr1.dat")
	vdat.printf(thr1,"%f\n")
	tdatfile.wopen("tdatfile.dat")
	tdat.printf(tdatfile,"%f\n")
    print spikes2.x[0]*dt-icdel //1st AP latency
	drwn1=1	
}else if(normcur==1.5 && drwn2==0){
    vdat.c.line(ana[num+1],tdat,color5,1)
	ana[num+1].flush
	thr1h.wopen("thr1h.dat")
	vdat.printf(thr1h,"%f\n")
 //   print spikes2.x[spiken-1]*dt-icdel 
	drwn2=1
	}else if(ic.amp==2){
    vdat.c.line(ana[num+1],tdat,color5,1)
	ana[num+1].flush
	thr2000.wopen("thr2000.dat")
	vdat.printf(thr2000,"%f\n")
//	print spikes2.x[spiken-1]*dt-icdel 
	drwn2=1
}else if(ic.amp==3){
   vdat.c.line(ana[num+1],tdat,color5,1)
	ana[num+1].flush
	thr3000.wopen("thr3000.dat")
	vdat.printf(thr3000,"%f\n")
//	print spikes2.x[spiken-1]*dt-icdel 
	drwn2=1
}else if(normcur==2 && drwn3==0){
    vdat.c.line(ana[num+2],tdat,color5,1)
	ana[num+2].flush
	thr2.wopen("thr2.dat")
	vdat.printf(thr2,"%f\n")
   if(spiken>1){
   print (spikes2.x[spiken-1]*dt-spikes2.x[spiken-2]*dt)/(spikes2.x[1]*dt-spikes2.x[0]*dt) // the last/1st of interspikeinterval(ISIlast/ISI1st)
   }else if(spiken==1){
  // print 1
   }
 //  print spikes2.x[spiken-1]*dt-icdel
	drwn3=1
}else if(normcur==2.5 && drwn4==0){
    vdat.c.line(ana[num+3],tdat,color5,1)
	ana[num+3].flush
	thr2h.wopen("thr2h.dat")
	vdat.printf(thr2h,"%f\n")
		drwn4=1
}else if(normcur==3 && drwn5==0){
    vdat.c.line(ana[num+4],tdat,color5,1)
	ana[num+4].flush
	thr3.wopen("thr3.dat")
	vdat.printf(thr3,"%f\n")
	drwn5=1
	
if(mode==1){
excute stoprun=1
}else if(mode==3){
		soma{
	gkcnqbar_kcnq1=gkcnqbar_kcnq1+100
	iint=threcur-10
	}
	if(gkcnqbar_kcnq1==2100){
	soma{
	gkcnqbar_kcnq1=0
	gkcnabar_kcna1=gkcnabar_kcna1+20
	iint=threcur0
	}
	}
	if(gkcnabar_kcna1==1020){
	excute stoprun=1
	}
	norm2=0
	run1()
}else if(mode==4){
soma{
	gkcnqbar_kcnq1=gkcnqbar_kcnq1+100
	iint=threcur-10
	}
	if(gkcnqbar_kcnq1==2100){
	soma{
	gkcnqbar_kcnq1=0
	gkcnab2bar_kcnab21=gkcnab2bar_kcnab21+20
	iint=threcur0
	}
	}
	if(gkcnab2bar_kcnab21==1020){
	excute stoprun=1
	}
	norm2=0
		run1()
}
color4=color4+1
color5=color5
}
}