strdef cmd, isomer
objref vbox1, vbox2, vbox3         // graphs variables
objref gv1, gi2, gi3               // graphs variables
objref volt_cl                     // variable for Point-Process
objref peak_vec 		   // vector for the peak current
objref v_vec			   // vector for the steps of voltage clamping

create soma
soma {
	diam = 50        // micron
	L    = 63.66198  // micron, so that area = 10000 micron2
	nseg = 1 	     // dimensionless
	cm   = 1         // uF/cm2
	Ra   = 70        // ohm-cm
	
	volt_cl = new VClamp_plus(.5)
}

y_min = -0.9 // variable for the maximum peak (graphic)
displ = 0    // current and voltage curves not showed (0), or 
		     // showed (1) during the simulation

incr = 5      		 // step increment
dens = 0	  		 // initializing clamping current
cl_st = -140
cl_end = 0

peak_vec = new Vector()
v_vec    = new Vector()

soma volt_cl=new VClamp_plus(.5)

proc Clamp() {

     peak_curr = 0     // inizializzazione del valore della corrente di picco
     volt_cl.amp[1] = $1			  // mV
     finitialize(hold_pot)
     
     while (t<tstop) {

             dens=volt_cl.i/area(.5)*100-soma.i_cap(.5) // clamping current in mA/cm2,
														// subtracted the capacitive current
             if (displ==1) {
				 gv1.line(t, soma.v(.5))
				 gv1.flush()
				 gi2.line(t, dens) 
				 gi2.flush()
             }
             
             if (t>(dur_st_cond+10)) {
             	
                if (abs(dens)>peak_curr) {
                	peak_curr=abs(dens)
            	}
             }
         fadvance()
     }
     
     peak_vec.append(peak_curr)
     v_vec.append($1)
     print peak_curr, $1   // debugging
     doEvents()
}

// erasing graphics 1 and 2

proc erase() {

     gv1.erase(0)
     gi2.erase(0)
     gv1.size(0, tstop, cl_st, cl_end)
     gi2.size(dur_st_cond+10, dur_st_cond+30, y_min, 0.1)
     gv1.beginline()
     gi2.beginline(2,1)

}

proc start() {
	
    peak_vec.resize(0)
    v_vec.resize(0)
    Tstop()

    for (i=cl_st; i<=cl_end; i=i+incr) {

		Clamp(i)

		erase()
	}
	
	peak_norm = peak_vec.max()
	for i=0, peak_vec.size()-1 {
		peak_vec.x[i]=peak_vec.x[i]/peak_norm
	}
	
    gi3.erase()
    gi3.size(cl_st-10, cl_end+10, 0, 1)
    gi3.begin()
	peak_vec.line(gi3,v_vec,1,1)
        
    for i=0, peak_vec.size()-1 {
		gi3.flush()
		doNotify()
    }
    
}

proc change_isomer() {
	num_to_strg(num_iso)
	sprint (cmd, "%s %s", "uninsert", isomer)
	execute(cmd)
	
	num_to_strg(nw_num_iso)
	sprint (cmd, "%s %s", "insert", isomer)
	execute(cmd)
	
	num_iso = nw_num_iso
}	
	
proc num_to_strg() {
	if ($1 == 1) {
		isomer = "na11a"
		y_min=-0.9
		hold_pot = -120      // holding potential
		dur_st_cond = 100    // duration of the conditioning stimulus
		dur_st_test = 20	 // test stimulus duration
		amp_st_test = -10	 // test stimulus amplitude
		cl_st = -140  		 // voltage-clamp starting value
		cl_end = 0    		 // voltage-clamp ending value
	}
	if ($1 == 2) {
		isomer = "na12a"
		y_min=-0.9
		hold_pot = -120      
		dur_st_cond = 100    
		dur_st_test = 20	 
		amp_st_test = -10	 
		cl_st = -140  		 
		cl_end = -10
	}
	if ($1 == 3) {
		isomer = "na13a"
		y_min=-0.9
		hold_pot = -90      
		dur_st_cond = 1000    
		dur_st_test = 20	 
		amp_st_test = -10	 
		cl_st = -100  		 
		cl_end = 15
	}
	if ($1 == 4) {
		isomer = "na14a"
		y_min=-0.8
		hold_pot = -120      
		dur_st_cond = 100    
		dur_st_test = 50	 
		amp_st_test = -10	 
		cl_st = -140  		 
		cl_end = -20
	}
	if ($1 == 5) {
		isomer = "na15a"
		y_min=-1.1
		hold_pot = -120      
		dur_st_cond = 500    
		dur_st_test = 20	 
		amp_st_test = -10	 
		cl_st = -120  		 
		cl_end = 0
	}
	if ($1 == 6) {
		isomer = "na16a"
		y_min=-1.7
		hold_pot = -90      
		dur_st_cond = 1000 
		dur_st_test = 20	 
		amp_st_test = 0	 
		cl_st = -120  		 
		cl_end = 0
	}
	if ($1 == 7) {
		isomer = "na17a"
		y_min=-1.5
		hold_pot = -140      
		dur_st_cond = 500    
		dur_st_test = 20	 
		amp_st_test = -20	 
		cl_st = -150  		 
		cl_end = -10
	}
	if ($1 == 8) {
		isomer = "na18a"
		y_min=-1.1
		hold_pot = -70      
		dur_st_cond = 500    
		dur_st_test = 40	 
		amp_st_test = 0	 
		cl_st = -80  		 
		cl_end = 20
	}
	if ($1 == 9) {
		isomer = "na19a"
		y_min=-3.3
		hold_pot = -120      
		dur_st_cond = 300    
		dur_st_test = 50	 
		amp_st_test = -40	 
		cl_st = -140  		 
		cl_end = 10
	}
}

proc Tstop() {
	 
	 volt_cl.dur[0] = 10	  		  // ms
     volt_cl.amp[0] = hold_pot   	  // mV
     volt_cl.dur[1] = dur_st_cond     // ms
     volt_cl.dur[2] = dur_st_test     // ms
     volt_cl.amp[2] = amp_st_test	  // mV
     volt_cl.dur[3] = 10		 	  // ms
     volt_cl.amp[3] = hold_pot		  // mV
     
     tstop = 0
     for i=0, 3 {
		tstop = tstop + volt_cl.dur[i]
	 }
}

access soma
num_iso = 1
nw_num_iso = num_iso
num_to_strg(num_iso)
sprint (cmd, "%s %s", "insert", isomer)
execute(cmd)
variable_domain("nw_num_iso", 1, 9)
	
{
xpanel("RunControl", 0)
hold_pot = -120
v_init = hold_pot
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 140
xvalue("t","t", 2 )
tstop = 140
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.0125
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 0.00999999
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(3,102)
}

{
xpanel("Start", 0)
xbutton("Start","start()")
celsius=22
xvalue("celsius")
ena=65
xvalue("ena")
xvalue("isomer (1-9)","nw_num_iso",1,"change_isomer()")
hold_pot=v_init
volt_cl.amp[0]=hold_pot
xvalue("holding_pot", "volt_cl.amp[0]")
volt_cl.dur[0]=10
xvalue("dur_pre_cond", "volt_cl.dur[0]", 1, "Tstop()")  // holding potential duration before conditioning stimulus
dur_st_cond = 100
xvalue("dur_st_cond", "dur_st_cond", 1, "Tstop()")   // conditioning stimulus duration
amp_st_test = -10	 
xvalue("amp_st_test", "amp_st_test")    // test stimulus amplitude
dur_st_test = 20
xvalue("dur_st_test", "dur_st_test", 1, "Tstop()")   // test stimulus duration
volt_cl.amp[3]=hold_pot
xvalue("amp_post_test", "volt_cl.amp[3]")
volt_cl.dur[3]=10
xvalue("dur_post_test", "volt_cl.dur[3]", 1, "Tstop()")  // holding potential duration after test stimulus
xvalue("starting_clamp", "cl_st")
xvalue("ending_clamp", "cl_end")
xvalue("incr", "incr")
xstatebutton("display", &displ)  // displaying single stimuli during simulation
xpanel(280,102)
}

// graphic for soma voltage

vbox1=new VBox()
vbox1.intercept(1)
gv1=new Graph()
gv1.size(0, tstop, cl_st, cl_end)
vbox1.intercept(0)
vbox1.map("Membrane voltage", 3, 600, -1, 0)

// graphic for clamp current (clamp.i)

vbox2=new VBox()
vbox2.intercept(1)
gi2=new Graph()
gi2.size(dur_st_cond+10, dur_st_cond+30, y_min, 0.1)
vbox2.intercept(0)
vbox2.map("Soma clamp current", 550, 10, 500, 350)

// grafico per corrente di inattivazione normalizzata

vbox3=new VBox()
vbox3.intercept(1)
gi3=new Graph()
gi3.size(cl_st-10, cl_end+10, 0, 1)
vbox3.intercept(0)
vbox3.map("Normalized current-voltage relation", 550, 400, 500, 350)