/*
DRG neuron from Kovalsky et at. : Simulation of Sensory Neuron Reveals a Key Role for Delayed Na+ Current
in Subthreshold Oscillation and Ectopic Discharge: Implications for Neuropathic Pain
J. Neurophysiology 102: 1430-1442, 2009

Contact: Marshall Devor, Ph.D. (marshlu@vms.huji.ac.il)
Department of Cell & Developmental Biology Institute of Life Sciences, 3-533  and Center for Research on Pain The Hebrew 
University of Jerusalem Jerusalem 91904, Israel
tel.: +972 2 6585085  FAX: +972 2 6586027

GUI control is added by Dongchul C. Lee, Ph.D. (Dongchul.Lee@bsci.com)
Principal Biomedical Engineer in Boston Scientific Neuromodulation, Valencia, CA,

*/

//Peak of IV_200 at -30 of 1.68 nA
//Osc start with stim=0.012, Vm=-57.21, freq=44.66, amp=2.3mV (fig 2A).
//stim=0.019 only two bursts (fig 2B). Vm=-57.14, osc freq=44.68.
//Range for the first isi in a burst in all the stimuli that generate bursts is 29.6-31.2 
//Wide bursts, stim=0.045 (fig 2C). Vm=-56.9, osc freq=45.55. Or stim=0.05, Vm=-56.8, osc freq=44.93.
//Spike height (small_stim=0.129) is 65.
//Tonic firing from stim=0.063, freq=36. Highest freq for tonic firing 99.3Hz (stim=1.362).
//Highest potential for osc is -36.63 (stim=1.62), osc freq = 109.77, amp=1.15.

load_file("nrngui.hoc")
wopen("out.txt")

steps_per_ms=100
dt=1/steps_per_ms
tstop=3000
secondorder=1
v_init=-56.8 // default: -57.42

create node
access node
node {  
	diam=50
	L=50		
	nseg=1			 
        cm=1

	insert Shh_6
	gnabar_Shh_6 = 25e-3 
	gl_Shh_6 = .00142
        el_Shh_6 = -65.5
	insert ls
        gnabar_ls = 1.28e-4
        insert inter
        gnabar_inter = 27e-3
        ena = 62   
        ek=-94
	insert kdr7
		gbar_kdr7 = 0
}

celsius=20

objectvar stim, small_stim
stim = new IClamp(0.5)
stim.del=100
stim.dur=2800
stim.amp=0.05 //default: 0.063

small_stim = new IClamp(0.5)
small_stim.del=100
small_stim.dur=1
small_stim.amp=0

objref vrec, inaVec, ilVec, ikVec, g ,g1, g2, g3

vrec = new Vector()
vrec.record(&v(0.5))
inaVec = new Vector()
ilVec = new Vector()
inaVec.record(&ina(0.5))
ilVec.record(&il_Shh_6(0.5))
ikVec = new Vector()
ikVec.record(&ik(0.5))

// GUI was added by Dongchul C. Lee
g = new Graph(0) // Voltage
g.addexpr("node.v(0.5)",1,1)
{g.view(0, -80, 3000, 100, 500, 0, 500, 150)}

g1 = new Graph(0) // M and H gate for fast Na+
g1.addexpr("node.m_Shh_6(0.5)",2,1)
g1.addexpr("node.h_Shh_6(0.5)",3,1)
{g1.view(0, 0, 3000, 1, 500, 250, 500, 150)}

g2 = new Graph(0) // M and H gate for intermediate Na+
g2.addexpr("node.m_inter(0.5)",2,1)
g2.addexpr("node.h_inter(0.5)",3,1)
{g2.view(0, 0, 3000, 1, 500, 500, 500, 150)}

g3 = new Graph(0) // M and H gate for slow Na+
g3.addexpr("node.m_ls(0.5)",2,1)
g3.addexpr("node.h_ls(0.5)",3,1)
{g3.view(0, 0, 3000, 1, 500, 750, 500, 150)}

graphList[0].append(g)
graphList[1].append(g1)
graphList[2].append(g2)
graphList[3].append(g3)

xopen("Main_Panels.ses")

/* Following lines were in original files.
run()	
inaVec.add(ilVec)
inaVec.add(ikVec)
for i=stim.del*steps_per_ms, (stim.del+stim.dur)*steps_per_ms {
	if(inaVec.x[i] < 0 && inaVec.x[i-1] >= 0) {
		//print i
	}
}
//computing osc freq:
first_peak=0
freq_sum=0
osc_num=0
sum_v_mid=0
osc_num2=0
for(i=170000; i <= 220000; i=i+1) {
	if(vrec.x[i] > vrec.x[i-1] && vrec.x[i] > vrec.x[i+1] && vrec.x[i] < -20) {
		if(first_peak == 0) {
			first_peak=i
		} else {
			osc_num=osc_num+1
			osc_freq = 1000/((i - first_peak)*dt)
			freq_sum=freq_sum+osc_freq
			print "osc_freq ", osc_freq
			first_peak=i
			//break
		}
	}	
	if(vrec.x[i] < vrec.x[i-1] && vrec.x[i] < vrec.x[i+1] && first_peak > 0) {
		//print "osc amp = ", vrec.x[first_peak] - vrec.x[i]
		v_mid = vrec.x[i] + (vrec.x[first_peak] - vrec.x[i])/2
		sum_v_mid = sum_v_mid + v_mid
		osc_num2 = osc_num2 + 1
		print v_mid
	}	
}
print freq_sum/osc_num
print sum_v_mid/osc_num2

for(j=0; j < 295000; j=j+10) {
	fprint("%f, %f\n", j, vrec.x[j])
}
wopen() */