// Firing simulator.hoc executes current-clamp recording of the model M-cell.
// Step or 500Hz pulse currents with 1pA or 20pA increments are given to the model M-cell.
//
// Takaki Watanabe
// wtakaki@m.u-tokyo.ac.jp

proc set_original() {
access soma
soma{
	insertoriginal()
	}
	gnabar_na = nstomho(gnabar_na1)
	gkhtbar_kht = nstomho(gkhtbar_kht1)
	gkabar_ka = nstomho(gkabar_ka1)
	g_leak = nstomho(g_leak1)
	gkcnqbar_kcnq = nstomho(gkcnqbar_kcnq1)
	gkcnabar_kcna = nstomho(gkcnabar_kcna1)
	gkcnab2bar_kcnab2 = nstomho(gkcnab2bar_kcnab21)
	vm0 = -85
}

proc iclamp2(){
icinvl1=0
icinvl=icinvl1
icdur=100
idelta=1
execute("objref ic, vc")
ic = new IClamp2(0.5)
ic.del=icdel
ic.dur=icdur
ic.invl=icinvl
modename="Step current"
mode=1
norm2=0
reGraph()
}

proc iclamp500(){
icinvl1=1.5
icinvl=icinvl1
icdur=0.5
idelta=20
execute("objref ic, vc")
ic = new IClamp500(0.5)
ic.del=icdel
ic.dur=icdur
ic.invl=icinvl
modename="500Hz pulse"
mode=2
norm2=0
reGraph()
}

proc run1(){
    access soma
    stoprun=0
	drwn1=0
	drwn2=0
	drwn3=0
	drwn4=0
	drwn5=0
	normcur=0
	celsius=25
	reGraph()
	ana[3].beginline(color5,1)
	ana[6].beginline(color5,1)
    ana[7].beginline(color5,1)
	ana[10].beginline(color5,1)		
	for i = 0,ncur-1 {
		icur.x[i] = iint+(i*(idelta))
		}		
	for i = 0, ncmd-1{
		vcmd.x[i] = vint+(i*(vdelta))
		}		
	tstop = 100
    gp[0].addexpr("soma.v(0.5)",1,1, 2*tstop,0,2)
    gp[1].addexpr("soma.ina(0.5)",1,1, 2*tstop,0,2)
    gp[1].addexpr("soma.ik(0.5)",2,1, 2*tstop,0,2)
    gp[2].addexpr("soma.gna_na(0.5)",1,1,0.8,0.9,2)
	gp[2].addexpr("soma.gkcna_kcna(0.5)",6,1,0.8,0.9,2)
	gp[2].addexpr("soma.gkcnab2_kcnab2(0.5)",6,1,0.8,0.9,2)
    gp[2].addexpr("soma.gkht_kht(0.5)",3,1,0.8,0.9,2)
	gp[2].addexpr("soma.gka_ka(0.5)",4,1,0.8,0.9,2) 
	gp[2].addexpr("soma.gkcnq_kcnq(0.5)",2,1,0.8,0.9,2)
    gp[3].addexpr("soma.i_cap(0.5)",2,1,0.8,0.9,2)
    gp[3].xexpr("soma.v(0.5)")
	for i=0,3{
	gp[i].begin()
	}

if (mode==1||mode==3||mode==4){
	for i = 0,ncur-1 {
	 if(normalize==0){
	ic.del=icdel
        ic.dur=icdur
        ic.invl=icinvl
		ic.amp = icur.x[i]/1000
	    v=vm0
		finitialize(v)
		fcurrent() 
		t=0
		while (t<tstop*1.5) {
		    fadvance()
				for j=0,3{
	gp[j].plot(t)
	}	
		}
		for j=0,3{
	gp[j].flush()
	}
		doNotify()
	
			for j=0,3{
	gp[j].begin()
	}
  run3()
	}else if (normalize==1){
	if(normcur==0){
	    ic.dur=icdur
        ic.invl=icinvl
		ic.amp = icur.x[i]/1000
	    v=vm0
		finitialize(v)
		fcurrent() 
		t=0
		while (t<tstop*1.5) {
		    fadvance()
				for j=0,3{
	gp[j].plot(t)
	}	
		}
		for j=0,3{
	gp[j].flush()
	}
		doNotify()
	
			for j=0,3{
	gp[j].begin()
	}
  run3()
  }else if(normcur > 0){
  for i=1,4{
	    ic.del=icdel
        ic.dur=icdur
        ic.invl=icinvl
		ic.amp = (1 + 0.5*i)*threcur/1000
	    v=vm0
		finitialize(v)
		fcurrent() 
		t=0
		while (t<tstop*1.5) {
		    fadvance()
				for j=0,3{
	gp[j].plot(t)
	}	
		}
		for j=0,3{
	gp[j].flush()
	}
		doNotify()
			for j=0,3{
	gp[j].begin()
	}
  normcur=1 + 0.5*i
  run3()

  }
  }
}
}
}else if(mode==2){
for i = 0,ncur-1 {
ic.del=icdel
        ic.dur=icdur
        ic.invl=icinvl
		ic.amp = icur.x[i]/1000
  if(ic.amp==3.02){
excute stoprun=1
}
	    v=vm0
		finitialize(v)
		fcurrent() 
		t=0
		while (t<tstop*1.5) {
		    fadvance()
				for j=0,3{
	gp[j].plot(t)
	}	
		}
		for j=0,3{
	gp[j].flush()
	}
		doNotify()
	
			for j=0,3{
	gp[j].begin()
	}
  run3()
}
}
}

proc delta20pA() {
	ncur =500
    iint = 0
    idelta=20
    icdel=5
	icinvl=icinvl1
    icamp=100
	xscale0=3
    xscale1=120
	mode=2
	normalize=1
	set_original()
	Clearrecord()
	}	
	
proc delta1pA()  {
	ncur =300
    idelta=1
    icdel=5
	icinvl=icinvl1
    icamp=100
	mode=1
	xscale0=3
    xscale1=120
	normalize=0
	}