OriginalDBS.hoc

//DBS of a GPi neuron (3D GPi dendrites, soma, axon)
//Matt Johnson (June 2007)
//NEURON 6.1 (Windows version)
//adapted from code by: Miocinovic 2006 and Gillies 2006

load_file("nrngui.hoc")
load_file("GPi_Fig4_cda_hocProps_Active.hoc")  //GPi cell geometry


//-----------------------------------------------------------------------------------------------------------
//  Load monkey-specific stimulation settings
//-----------------------------------------------------------------------------------------------------------

monkey = 3  //momo = 1; mary = 2; nela = 3; raisin=4;

strdef waveform_filename  //stimulus waveform made with capacitive Femlab model and fourier_waveform.m
strdef voltage_filename   //voltage file (voltages for each cell and each compartment) made with dbs_pop_coord2.m
strdef output_file1, output_file2

if (monkey == 1) { 
	
	//MOMO (bipo3c1a 2V 3.5V 5.5V) & (bipo0c2a 2V 3.5V 5.5V)
	waveform_filename = "ffw_reconstr_bipo3c1a_encap250_136Hz_90us.txt" //"waveform_bargad_mono_biphasic_200us_136Hz.txt"
	voltage_filename = "NEURON_0d43V_momo_gpi_bipo3c1a_pop1000.txt"     //"GPi_Fig4new_mono1c_pop63_1V.txt"
	output_file1 = "Momo_GPi_bipo3c1a_pop000-099_threshV.dat" 
	output_file2 = "Momo_GPi_bipo3c1a_pop000-099_threshV.log"

	//Extracellular stimulation parameters
	Vamp = 2 	//[V] extracell. stimulation voltage
	pw = 0.09  	//ms
	
}

if (monkey == 2) {
	
	//MARY bipo 2V 3V 3.5V
	waveform_filename = "ffw_reconstr_bipo3c1a_encap250_136Hz_90us.txt"
	voltage_filename = "NEURON_0d43V_mary_gpi_bipo1c0a_pop1000.txt"
	output_file1 = "Mary_GPi_bipo1c0a_pop000-099_6V.dat" 
	output_file2 = "Mary_GPi_bipo1c0a_pop000-099_6V.log"

	//Extracellular stimulation parameters
	Vamp = 13.95  //[V] extracell. stimulation voltage
	pw = 0.09    //ms  //LOAD CORRECT FOURIER-DERIVED WAVEFORM 

}

if (monkey == 3) {
	
	//NELA mono 2V
	waveform_filename = "ffw_reconstr_bipo3c1a_encap250_136Hz_90us.txt"
	voltage_filename = "NEURON_0d43V_nela_gpi_bipo3c1a_pop1000.txt"
	output_file1 = "Nela_GPi_bipo3c1a_pop000-099_3V.dat" 
	output_file2 = "Nela_GPi_bipo3c1a_pop000-099_3V.log"

	//Extracellular stimulation parameters
	Vamp = 3  //[V] extracell. stimulation voltage
	pw = 0.09    //ms  //LOAD CORRECT FOURIER-DERIVED WAVEFORM 

}

if (monkey == 4) {
	
	//RAISIN mono 2V
	waveform_filename = "ffw_reconstr_bipo3c1a_encap250_136Hz_90us.txt"
	voltage_filename = "NEURON_0d43V_raisin_gpi_mono2c_pop1000.txt"
	output_file1 = "Raisin_GPi_mono2c_pop000-099_threshV.dat" 
	output_file2 = "Raisin_GPi_mono2c_pop000-099_threshV.log"

	//Extracellular stimulation parameters
	Vamp = 2  //[V] extracell. stimulation voltage
	pw = 0.09    //ms  //LOAD CORRECT FOURIER-DERIVED WAVEFORM 

}

if (monkey == 5) {
	
	//Model for JNP paper
	waveform_filename = "ffw_reconstr_bipo3c1a_encap250_136Hz_90us.txt"
	voltage_filename = "GPi_Fig4new_mono1c_pop63_1V.txt"
	output_file1 = "GPi_Fig4_mono1c_2V_temp.dat" 
	output_file2 = "GPi_Fig4_mono1c_2V_temp.log"

	//Extracellular stimulation parameters
	Vamp = 0.5  //[V] extracell. stimulation voltage
	pw = 0.09    //ms  //LOAD CORRECT FOURIER-DERIVED WAVEFORM 

}

//-----------------------------------------------------------------------------------------------------------
//  Set up stimulation and simulation parameters
//-----------------------------------------------------------------------------------------------------------

//Extracellular stimulation parameters
delay = 50		//ms  was 50 (145 for synapse testing)
freq = 136 		//Hz
polarity = -1 	//-1 for cathodic, +1 for anodic , -1 for bipolar
ratio = 10     	//pseudomonophasic: amplitude amp/ratio, duration pw*ratio (and polarity -1*polarity)
num_test_pulses = 10
dur = num_test_pulses*1000/freq  //ms

//Simulation parameters
dt = 0.01               //ms
tstop = delay+dur+50  	//ms
steps_per_ms = 100

//Voltage parameters
max_pop = 1000 				//max number of cells in population -- was 60
objref V_raw, Ve, V
V_raw = new Vector(total*max_pop,0)	//total number of compartments (comes from geometry file) 
Ve = new Vector(total,0)
V = new Vector(total,0)

//Synapse input vector trains
objref VAMPAdend, VGABAadend, VAMPAsoma, VGABAasoma
objref VAMPAdende, VGABAadende, VAMPAsomae, VGABAasomae
VAMPAdend = new Vector(numAMPA_dend,0)
VGABAadend = new Vector(numGABAa_dend,0)
VAMPAsoma = new Vector(numAMPA_soma,0)
VGABAasoma = new Vector(numGABAa_soma,0)
VAMPAdende = new Vector(numAMPA_dend,0)
VGABAadende = new Vector(numGABAa_dend,0)
VAMPAsomae = new Vector(numAMPA_soma,0)
VGABAasomae = new Vector(numGABAa_soma,0)

//Variables for calculating average firing activity at the end of the axon
avg_inst_freq = 0
std_inst_freq = 0
avg_period = 0
std_period = 0
avg_freq = 0
objref order_freq
order_freq = new Vector (dur, 0)   //instantaneous firing rate for each spike
dd = 0 //counter for calculating freq info (number of interspike intervals)



//-----------------------------------------------------------------------------------------------------------
//  Set up xyz scale bars
//-----------------------------------------------------------------------------------------------------------

create xScale, yScale, zScale
proc anatscale() {
	if ($4>0) {  // if length arg is <= 0 then do nothing
		xScale {
			pt3dclear()
			pt3dadd($1, $2, $3, 1)
			pt3dadd($1+$4, $2, $3, 1)
		}
		yScale {
			pt3dclear()
			pt3dadd($1, $2, $3, 1)
			pt3dadd($1, $2+$4, $3, 1)
		}
		zScale {
			pt3dclear()
			pt3dadd($1, $2, $3, 1)
			pt3dadd($1, $2, $3+$4, 1)
		}
	}
}
anatscale(0,0,0,100)  // origin xyz and length



//-----------------------------------------------------------------------------------------------------------
// Create an extracellular stimulus clamp
//-----------------------------------------------------------------------------------------------------------

objref exIClmp
exIClmp = new IClamp()
proc stimulus(){
	
	//extracellular stimulus
	Vstim = Vamp * polarity
	electrode {                 //electrode is created in morphology file 
        exIClmp.loc(0.5)
		exIClmp.del = 0         //was 0, onset delay, will use vectorplay 
    	exIClmp.dur = 500000  	//was 1e6, duration, will use vectorplay 
		exIClmp.amp = 0         //was 0, will use vectorplay to determine stimulus amplitude at any time point
	}

	for i=0,total-1 {
		Ve.x[i]=Vstim*V.x[i]*1e3	//mV
	}


    for i=0,numAMPA_dend-1 {
        VAMPAdende.x[i]=Vstim*VAMPAdend.x[i]*3e2
    }
    for i=0,numGABAa_dend-1 {
        VGABAadende.x[i]=Vstim*VGABAadend.x[i]*3e2
    }
    for i=0,numAMPA_soma-1 {
        VAMPAsomae.x[i]=Vstim*VAMPAsoma.x[i]*3e2
    }
    for i=0,numGABAa_soma-1 {
        VGABAasomae.x[i]=Vstim*VGABAasoma.x[i]*3e2
    }

}
stimulus()



//-----------------------------------------------------------------------------------------------------------
// Define a custom (fourier-derived) waveform for the extracellular stimulus
//-----------------------------------------------------------------------------------------------------------

objref extra_stim
time_step = int(tstop/dt) + 1   
extra_stim = new Vector(time_step,-1)  //RESIZE IF CHANGING TSTOP LATER

//Load one cycle of fourier-derived waveform shape (normalized to 1V; dt=0.01); loads vector fdw
objref fdw
cycle_size = int(1000/(freq*dt)) + 1
fdw = new Vector(cycle_size,0)   
xopen(waveform_filename)

//Calculate extra_stim values for each time-step and feed into exIClmp.amp
proc calc_extra_stim(){ local tmp_del, tmp_dur, my_time, i, j

    my_time = 0
	if (freq == 0 || dur == 0 || Vamp == 0) {                   //no stimulation
    	for i = 0, time_step-1 { 
           	extra_stim.x[i] = 0
        }
    } else {
        i = 0
        while(my_time < delay && i < extra_stim.size()) {  		//before stimulus train
           	extra_stim.x[i] = 0
           	i = i + 1
	      	my_time = my_time + dt      
        }
        while (my_time < delay+dur && i < extra_stim.size()) { 	//during stimulus train
            for j = 0, cycle_size-1 {
             	extra_stim.x[i] = fdw.x[j]   //extra_stim has no units
             	i = i + 1
                my_time = my_time + dt 
            }        
        }
	  	while (i < extra_stim.size()) {     				//after stimulus train  
            extra_stim.x[i] = 0
            i = i + 1
            my_time = my_time + dt         
      	}	
    }
}
calc_extra_stim()



//-----------------------------------------------------------------------------------------------------------
// Create an action potential counter
//-----------------------------------------------------------------------------------------------------------

objref apc, apc_times, apc_soma, apc_times_soma
proc setup_AP_count(){

	AP_ct_node = -1 
	node[26] apc = new APCount(0.2)
	AP_ct_node = 26
	apc_times = new Vector()
	apc.thresh = -20 //mV
	apc.record(apc_times)

	soma[3] apc_soma = new APCount(0.2)   
	apc_times_soma = new Vector()
	apc_soma.thresh = -20 //mV
	apc_soma.record(apc_times_soma)

}
setup_AP_count()

//Set up AP counters in every node and branches (if axon has any)
num_ap_counters = axonnodes+1

//Set up AP counters in every node
objectvar AP_counters[num_ap_counters]  //AP counter at every node
objectvar AP_timerec[num_ap_counters]   //records AP time at every node
for i = 0, num_ap_counters-1 {
    AP_timerec[i] = new Vector()
}

proc setup_APc_all(){
	for i = 0,axonnodes-1 {
		node[i] AP_counters[i] = new APCount(.5)    //put AP counter in every node
		AP_counters[i].record(AP_timerec[i])        //records times of APs
	}
	soma[1] AP_counters[num_ap_counters-1] = new APCount(.5)  //put AP counter in soma //MULTICOMPARTMENT SOMA
	AP_counters[num_ap_counters-1].record(AP_timerec[num_ap_counters-1])      

}
setup_APc_all()



//-----------------------------------------------------------------------------------------------------------
// Determines if AP occurred
//-----------------------------------------------------------------------------------------------------------

//Threshold seeking functions
objref extra_APs, required_APs
required_APs = new Vector(1000,0) 	//times of current stimulus pulses
num_req_APs = 0                     //number of APs expected in response to stimulus pulses
extra_APs = new Vector(1000,0)      //times of APs during stimulation not classified as artifact
num_extra = 0                       //size of extra_APs vector
cyc = 0                             //number of stimulation pulses
epsilon = 0.01                      //was 0.1 //was 0.001 //allowed error when looking for threshold
AP_delay = 3                        //was 4 ms, time allowed between stimulus pulse initiation and AP at the 
                                    //AP_counter to prevent considering AP as part of the stimulus
num_req = 0                         //number of required APs that happened 
soma_orig_ap = 0                    //number of APs that originated in soma but counted as stimulus related
soma_orig_ap2 = 0                   //number of APs that originated in soma but counted as extra APs


//Calculate times of current stimulus pulses (starting time), and store them in Vector required_APs
//APs should follow after a certain time whose maximum is AP_delay
proc calculate_required_APs() {  local i

	num_req_APs = 0

    if (freq == 1) {
		cyc = 1
    } else {
       	if (dur > tstop) {
           	cyc = int((tstop-delay)*freq/1000) 
        } else {
            cyc = int(dur*freq/1000) 
        }
    }

    for i = 0, cyc-1 {
        required_APs.x[num_req_APs] = delay+i*1000/freq
        num_req_APs = num_req_APs + 1
    } 

    if (num_req_APs != num_test_pulses) {
        print "Freq = ", freq, "...Number of required APs is not ", num_test_pulses, "!!, but ", num_req_APs
    }
}

//Checks if cell responds to all stimulus pulses and records any extra APs to a vector extra_APs
func response_to_stimulus(){  local i, j, flag, flag2

  	num_req = 0 
    num_extra = 0
    flag = 0     
    flag2 = 0
    soma_orig_ap = 0
    soma_orig_ap2 = 0

    for i = 0, apc.n-1 {
        if ((apc_times.x[i] > delay) && (apc_times.x[i] <= delay+dur+AP_delay)) { //AP_delay added if there are any APs that are responding to stimulus pulse at the very end of pulse duration
            flag = 0  
            for j = flag2, num_req_APs-1 {
                if ((apc_times.x[i] > required_APs.x[j]) && (apc_times.x[i] <= required_APs.x[j]+AP_delay)) {
                    if (soma_ap(i)==1){
                        soma_orig_ap = soma_orig_ap + 1
                    }
                    num_req = num_req + 1
                    flag = 1
                    flag2 = j + 1
                    break
                }  
            }
            if (flag == 0) {
		        if (soma_ap(i)==1){
                    soma_orig_ap2 = soma_orig_ap2 + 1
                }
                extra_APs.x[num_extra] = apc_times.x[i]
                num_extra = num_extra + 1
            }
        }  
    }

   	print "      Num required current pulses = ", num_req_APs, "  Num required that happened = ", num_req
   	print "      Num_extra = ", num_extra

    // if (num_req > 0) {  //there is at least one AP in response to stimulation
    // if (num_req == num_req_APs) {  //there is AP following every stimulus pulse
	if (num_req >= 0.8*num_req_APs) {  //there is AP following 80% of stimulus pulses
      	return 1
    }

    return 0
}          

//returns 1 is AP originated in soma, 0 otherwise 
//ap_number is ordinal number of AP
func soma_ap(){local soma_ap, node0_ap

	ap_number=$1
    if ((AP_counters[num_ap_counters-1].n) <= ap_number || (AP_counters[0].n) <= ap_number) {
        return 0 //couldn't determine if AP originated in soma
    }

    soma_ap = AP_timerec[num_ap_counters-1].x[ap_number]
	node0_ap = AP_timerec[0].x[ap_number]

    //print "Soma ap time = ", soma_ap, "Node 0 AP time = ", node0_ap

    if (soma_ap < node0_ap) {
        return 1
    }

    return 0
}


//Run trial
func trial() {
	Vamp = $1   //in volts
    print "                 Trial voltage = ", Vamp*polarity, " V"
    stimulus()
	init() 
	run()
    if (response_to_stimulus() == 1) {
        return 1
    }
    return 0
}

//Return voltage threshold in V
func threshold() {
	strength = 6   //was 0.5 in Volts
    low_limit = 1e-3
	lbound = low_limit
    up_limit = 7    // was 50
	ubound = up_limit
	while (lbound == low_limit || ubound == up_limit) {
		excited = trial(strength)
		if (excited > 0) {
			ubound = strength
			strength = ubound/2
		} else {
            lbound = strength
		    strength = 2*lbound
        }
        if (lbound>ubound)  return up_limit  
        if (strength>ubound)  return up_limit
	}
	strength = (ubound+lbound)/2
    while((abs(ubound-lbound))>(abs(epsilon*ubound))) {
		excited = trial(strength)
		if (excited > 0) {
			ubound = strength
			strength = (ubound+lbound)/2
		} else {
			lbound = strength
			strength = (3*ubound+lbound)/4
		}
	}
    //return strength  // threshold in Volts
    //strength is between ubound and lbound; ubound is returned because that is the 
    //closest value to strength for which we know causes an AP (strength might not be enough)  
    return ubound   
}

//-----------------------------------------------------------------------------------------------------------
// Find AP initiation site
//-----------------------------------------------------------------------------------------------------------

objref temp_vec
temp_vec = new Vector(num_ap_counters,0)
proc AP_init_site(){

    AP_site1 = -1   	//node where AP appears first
    AP_site2 = -1   	//node where AP appears second
	time1 = -1          //time of AP appearance at site1
	time2 = -1          //time of AP appearance at site2

    for qw = 0, num_ap_counters-1 {
        temp_vec.x[qw] = 100000  //set it to a big number so it does not interfere with finding shortest time
        for wh = 0, AP_counters[qw].n-1 {
            //delay+pw*ratio  first AP when stim starts, but after first prepulse
			if (AP_timerec[qw].x[wh] > delay) {   //first AP when stim starts
                temp_vec.x[qw] = AP_timerec[qw].x[wh]
				break
			}
		}
    } 

    if (temp_vec.min() ==  100000) {
        print "NO AP during stimulation...."
    } else {
        //find index of AP counter that recorded shortest time (ie where AP first appeared)
        AP_site1 = temp_vec.min_ind()
        time1 = temp_vec.x[AP_site1]

	    //find second site for AP initialization
        temp_vec.x[AP_site1]= 100000
        AP_site2 = temp_vec.min_ind()  //2nd most depolarized node at time1
        time2 = temp_vec.x[AP_site2]
    }
}

//-----------------------------------------------------------------------------------------------------------
// Calculate the firing rate, average instantaneous rate, and firing period
//-----------------------------------------------------------------------------------------------------------

//Consider time between time_beg and time_end (given as arguments); AP counter is distal node
proc calculate_freq(){ local time_beg, time_end, sum1, sum2, sum3, sum4, dd, kk, gg

    time_beg = $1
    time_end = $2

    sum1 = 0
    sum2 = 0
    sum3 = 0
    sum4 = 0
    dd = 0  	//number of interspike intervals
    gg = 0  	//number of spikes

   	for kk = 0, apc.n-2 {
        if ((apc_times.x[kk] > time_beg) && (apc_times.x[kk+1] <= time_end)) {
            //    print "1st = ", apc_times.x[kk], "  2nd = ", apc_times.x[kk+1]
            tmp5 = apc_times.x[kk+1]-apc_times.x[kk]
            sum1 = sum1 + tmp5
            sum2 = sum2 + tmp5^2
            sum3 = sum3 + (1000/tmp5)   //convert from ms to seconds, then to Hz
            sum4 = sum4 + (1000/tmp5)^2
            order_freq.x[dd] = (1000/tmp5)   //convert from ms to seconds, then to Hz
            dd=dd+1
        }
    }

    avg_period = 0
    avg_inst_freq = 0
    if (dd != 0)  {
        avg_period = sum1/dd
        avg_inst_freq = sum3/dd
    }
    for kk = 0, apc.n-1 {
        if ((apc_times.x[kk] > time_beg) && (apc_times.x[kk] <= time_end)) gg = gg+1
    }	 
    avg_freq = 1000*gg/(time_end-time_beg)  //convert time from ms to sec

    std_period = 0 
    std_inst_freq = 0
    if (dd-1 > 0)  {
        var_period = (sum2 - sum1^2/dd)/(dd-1)
        var_freq = (sum4 - sum3^2/dd)/(dd-1)
        if (var_period > 0) {
            std_period = sqrt(var_period) 
        } 
        if (var_freq > 0) {
            std_inst_freq = sqrt(var_freq) 
        }          
    }
}

//calculate firing freq (#spikes/recording_time), average instantaneous firing freq, firing period (and standard deviations)
//consider time between time_beg and time_end (given as arguments)
//use AP counter at the soma 
proc calculate_freq_soma(){ local time_beg, time_end, sum1, sum2, sum3, sum4, dd, kk, gg

    time_beg = $1
    time_end = $2

    sum1 = 0
    sum2 = 0
    sum3 = 0
    sum4 = 0
    dd = 0  //number of interspike intervals
    gg = 0 //number of spikes
    for kk = 0, apc_soma.n-2 {
        if ((apc_times_soma.x[kk] > time_beg) && (apc_times_soma.x[kk+1] <= time_end)) {
            //    print "1st = ", apc_times_soma.x[kk], "  2nd = ", apc_times_soma.x[kk+1]
            tmp5 = apc_times_soma.x[kk+1]-apc_times_soma.x[kk]
            sum1 = sum1 + tmp5
            sum2 = sum2 + tmp5^2
            sum3 = sum3 + (1000/tmp5)   //convert from ms to seconds, then to Hz
            sum4 = sum4 + (1000/tmp5)^2
             
            order_freq.x[dd] = (1000/tmp5)   //convert from ms to seconds, then to Hz
            dd=dd+1
        }
    }

    avg_period_soma = 0
    avg_inst_freq_soma = 0
    if (dd != 0)  {
        avg_period_soma = sum1/dd
        avg_inst_freq_soma = sum3/dd
    }
    for kk = 0, apc_soma.n-1 {
        if ((apc_times_soma.x[kk] > time_beg) && (apc_times_soma.x[kk] <= time_end)) gg = gg+1
    }	 
    avg_freq_soma = 1000*gg/(time_end-time_beg)  //convert time from ms to sec
 
    std_period_soma = 0 
    std_inst_freq_soma = 0
    if (dd-1 > 0)  {
        var_period = (sum2 - sum1^2/dd)/(dd-1)
        var_freq = (sum4 - sum3^2/dd)/(dd-1)
        if (var_period > 0) {
            std_period_soma = sqrt(var_period) 
        } 
        if (var_freq > 0) {
            std_inst_freq_soma = sqrt(var_freq) 
        }          
    }
}

//Calculate time between end of stimulus pulse and first spontaneous spike
//Return -1 if there are no spikes after stimulus ended
//Wait wait_time milliseconds after stimulus ends before looking for an after stimulus spike
func calculate_delay(){ local wait_time
    
    wait_time = 2 	//ms   //WAS 2

    for kk = 0, apc.n-1 {
        if (apc_times.x[kk] > delay+dur+wait_time) {
            return apc_times.x[kk]-(delay+dur)
        }
    }

    return -1
}

//Calculate instantaneous firing frequency of first interspike interval during stimulus 
//Return -1 if there are no spikes (or only one) during stimulus
func calculate_first_int(){ 
    
	for kk = 0, apc.n-2 {
     	if ((apc_times.x[kk] > delay) && (apc_times.x[kk+1] <= delay+dur)) {
           	return 1000/(apc_times.x[kk+1]-apc_times.x[kk])  //convert to seconds, then to Hz
      	}
    }
   
    return -1
}

//Procedure that returns number of APs between time time_beg and time_end
func num_APs() { local count 

	time_beg = $1
  	time_end = $2

	count = 0
	for kk = 0, apc.n-1 {
		if ((apc_times.x[kk] > time_beg) && (apc_times.x[kk] <= time_end)) {
	 		count = count + 1
		}
	}
	return count
}


//-----------------------------------------------------------------------------------------------------------
// Initialize the simulation
//-----------------------------------------------------------------------------------------------------------

proc advance() {

	for i=0,total-1 {
		s[i].sec.e_extracellular(0.5)=(exIClmp.i)*Ve.x[i]	//mV//
	}

	for i=0,numAMPA_dend-1 {
		affAMPAdend[i].sec.e_extracellular(0.5)=(exIClmp.i)*VAMPAdende.x[i]	//mV//
	}
	for i=0,numGABAa_dend-1 {
		affGABAadend[i].sec.e_extracellular(0.5)=(exIClmp.i)*VGABAadende.x[i]	//mV//
	}
	for i=0,numAMPA_soma-1 {
		affAMPAsoma[i].sec.e_extracellular(0.5)=(exIClmp.i)*VAMPAsomae.x[i]     //mV//
	}
	for i=0,numGABAa_soma-1 {
		affGABAasoma[i].sec.e_extracellular(0.5)=(exIClmp.i)*VGABAasomae.x[i]	//mV//
	}

	//use to stop threshold seeking if after 20% of stimulus pulses there aren't at least 10% of required APs
	AP_during_stim = num_APs(delay,t)
    //if ( (t > delay+(1000*(0.20*num_test_pulses)/freq)+AP_delay) && (AP_during_stim < 0.10*num_req_APs)) {
	//	print "STOP RUN at t=", t, " because there were ", AP_during_stim, " APs after ", 0.20*num_test_pulses," test pulses"
	//	stoprun=1
	//}

	fadvance()
}

extra_stim.play(&exIClmp.amp, dt)

finitialize(v_init)
fcurrent()
xopen("GPi_dbs_fem.ses")    
//xopen("GPi_Fig4_synapseprops.ses")    

//-----------------------------------------------------------------------------------------------------------
// Load the extracellular voltage field data
//-----------------------------------------------------------------------------------------------------------

objref cell_pos
cell_pos=new Vector(max_pop*3,0) 

xopen(voltage_filename)  //load extracell. voltages from FEMLAB

//locations of neurons 
//num_cells is number of cells in population (loaded with voltage data)
objectvar cell_coords[num_cells]   //coordinates of cells in population (their center coords)
for i = 0, num_cells-1 {
    cell_coords[i] = new Vector(3,0)
    cell_coords[i].x[0] = cell_pos.x[i*3]   //cell_pos is also loaded with voltage data
    cell_coords[i].x[1] = cell_pos.x[i*3+1]
    cell_coords[i].x[2] = cell_pos.x[i*3+2]
}

//Extracellular stimulus pulse frequency (Hz)
objref freq_vec 
freq_vec = new Vector (1,0)
freq_vec.x[0]=136

//-----------------------------------------------------------------------------------------------------------
// Start the simulation
//-----------------------------------------------------------------------------------------------------------

print "Simulation running...\n"

objref f1, f2
f1 = new File()
f2 = new File()

f1.wopen(output_file1)  
f2.wopen(output_file2)

f1.printf("DBS voltage pulse train stimulation of axon POPULATION\n\n")
f2.printf("DBS voltage pulse train stimulation of axon POPULATION\n\n")

f1.printf("Stimulus waveform used (fourier-derived): %s\n", waveform_filename)
f1.printf("Voltage file used : %s\n", voltage_filename)
f1.printf("Voltage stimulus params: freq=variable Hz, pw=%f, delay=%f ms, dur=variable, amp=%f, polarity=%d, # test pulses = %d \n", pw, delay, Vamp, polarity, num_test_pulses)
f1.printf("Other params: tstop=%f ms, dt=%f ms, steps_per_ms=%f ms, AP_delay=%f, number of nodes=%d, epsilon=%f, AP_counter at node %d, ratio=%f\n", tstop, dt, steps_per_ms, AP_delay, axonnodes, epsilon, AP_ct_node, ratio)

f2.printf("Voltage stimulus params: freq=variable Hz, pw=%f, delay=%f ms, dur=variable, amp=%f, polarity=%d, # test pulses = %d \n", pw, delay, Vamp, polarity, num_test_pulses)
f2.printf("Other params: tstop=%f ms, dt=%f ms, steps_per_ms=%f ms, AP_delay=%f, number of nodes=%d, epsilon=%f, AP_counter at node %d \n", tstop, dt, steps_per_ms, AP_delay, axonnodes, epsilon, AP_ct_node)

f2.printf("Cell locations (x,y,z)\n")
for i = 0, num_cells-1 {
    f2.printf("%f\t%f\t%f\n", cell_coords[i].x[0], cell_coords[i].x[1], cell_coords[i].x[2])
}

if (dt != 0.01) {
	print "dt is not 0.01!"
}

//Main loop of the program
for ff = 0,0 {  //0 to 4

    print "Stim_freq = ", freq
    f1.printf("\nStimulus voltage pulse frequency: %f\n", freq)
    f2.printf("\n\nStimulus voltage pulse frequency: %f\n", freq)

    calculate_required_APs()

    f2.printf("Times of extracell. stimulus pulses\n")
    f2.printf("Voltage Pulse \t Time \n")
    for i = 0, num_req_APs-1 {
        f2.printf("%d \t %f\n", i+1, required_APs.x[i])   //print times of voltage extracell. stimulus pulses
    } 

    f1.printf("Cell # \t X \t Y \t Z \t Thresh voltage(V) \t ubound \t # required APs \t #extra APs during stim \t")
    f1.printf("Avg spikes/sec B(Hz) \t Avg Inst firing B (Hz) \t Inst rate STD B (Hz) \t Firing period B (ms) \t Firing period STD B (ms) \t ")
    f1.printf("Avg spikes/sec D (Hz) \t Avg Inst firing D \t Inst rate STD D (Hz) \t Firing period D (ms) \t Firing period STD D (ms) \t ")
    f1.printf("Avg spikes/sec A (Hz) \t Avg Inst firing A \t Inst rate STD A (Hz) \t Firing period A (ms) \t Firing period STD A (ms) \t ")
    f1.printf("Avg freq D aft 100ms(Hz) \t Delay after stim (ms) \t")
    f1.printf("AP site#1 \t AP site#2 \t AP time1 \t AP time2\n")

    //go thru different cells
    for jj = 0,99 {      //num_cells-1 {

        print "     Cell# ", jj, ", Cell location ", cell_coords[jj].x[0], cell_coords[jj].x[1], cell_coords[jj].x[2]
        f2.printf("\n\nCell position: %d %d %d\n", cell_coords[jj].x[0], cell_coords[jj].x[1], cell_coords[jj].x[2])

        for kk = 0, total-1 {
            V.x[kk] = V_raw.x[jj*total+kk]
        }

        for kk = 0, numAMPA_dend-1 {
            VAMPAdend.x[kk] = V_raw.x[jj*total+kk]
        }
        for kk = 0, numGABAa_dend-1 {
            VGABAadend.x[kk] = V_raw.x[jj*total+kk]
        }
        for kk = 0, numAMPA_soma-1 {
            VAMPAsoma.x[kk] = V_raw.x[jj*total+kk]
        }
        for kk = 0, numGABAa_soma-1 {
            VGABAasoma.x[kk] = V_raw.x[jj*total+kk]
        }

        //ubound = threshold()   //also UNCOMMENT portion of advance()
        //thresh_curr = ubound
        thresh_curr = -20
        ubound = Vamp //in Volts
        trial(ubound)  //do it again to get AP times for determined threshold 

        calculate_freq(delay, delay+dur)  //during stimulation
        calculate_freq_soma(delay, delay+dur)  //during stimulation
        AP_init_site()

        print " During stim --> avg_fr =" , avg_freq, ", avg_inst_fr =" , avg_inst_freq, "std_inst_fr = ", std_inst_freq 
        f1.printf("%d\t %f\t %f\t %f\t %f\t %f\t %d\t", jj,cell_coords[jj].x[0], cell_coords[jj].x[1], cell_coords[jj].x[2],thresh_curr, ubound, num_req)
        f1.printf("%d\t %f\t %f\t %f\t %f\t %f\t %d\t %d\t %f\t %f\t %d\t %d\t", num_extra, avg_freq, avg_inst_freq, std_inst_freq, avg_freq_soma, avg_inst_freq_soma, AP_site1, AP_site2, time1, time2,soma_orig_ap,soma_orig_ap2 )

        //f1.printf("%d \t %f \t %f \t %f \t %f \t %f \t %d\t %d\t", jj,cell_coords[jj].x[0], cell_coords[jj].x[1], cell_coords[jj].x[2],thresh_curr, ubound, num_req, num_extra)

        //before stimulation
        calculate_freq(0, delay)
        print "Before stimulation"
        print "     avg_fr =" , avg_freq, "Avg Inst fr = ", avg_inst_freq, "Hz, std_fr = ", std_inst_freq
        f1.printf("%f \t %f \t %f \t %f \t %f \t ", avg_freq, avg_inst_freq, std_inst_freq, avg_period, std_period)

        //during stimulation
        calculate_freq(delay, delay+dur)
        print "During stimulation"
        print "     avg_fr =" , avg_freq, "Avg Inst fr = ", avg_inst_freq, "Hz, std_fr = ", std_inst_freq
        f1.printf("%f \t %f \t %f \t %f \t %f \t ", avg_freq, avg_inst_freq, std_inst_freq, avg_period, std_period)

        //after stimulation
        calculate_freq(delay+dur, tstop)
        print "After stimulation"
        print "     avg_fr =" , avg_freq, "Avg Inst fr = ", avg_inst_freq, "Hz, std_fr = ", std_inst_freq
        f1.printf("%f \t %f \t %f \t %f \t %f \t ", avg_freq, avg_inst_freq, std_inst_freq, avg_period, std_period)
        spike_del = calculate_delay()
        print "     spike delay =" , spike_del, " ms"

        calculate_freq(delay+100, delay+dur)  //100ms after stim begins
        f1.printf("%f \t %f \t ", avg_freq, spike_del)
        f1.printf("%d \t %d \t %f \t %f\n", AP_site1, AP_site2, time1, time2)

        f2.printf("\n\nAll AP times for cell (at axon end)# %d(ms)\n", jj)
        for kk=0, apc.n-1 {
            f2.printf("\t %f\n", apc_times.x[kk])
        }
        f2.printf("Extra AP times (at end of axon) during stimulation (ms) of cell# %d\n", jj)
        for kk=0, num_extra-1 {
            f2.printf("\t %f\n", extra_APs.x[kk])
        }
        f2.printf("Soma AP times all (ms) of cell# %d\n", jj)
        for kk=0, AP_counters[num_ap_counters-1].n-1 {
            f2.printf("\t %f\n", AP_timerec[num_ap_counters-1].x[kk])
        }	

        print "jj = ", jj

    }
    print "jj = ", jj
    print "num_cells = ", num_cells
}

f1.close()
f2.close()