// Created by Brandon Thio 7 February 2018
begintemplate cFiberBuilder
// builds different c-fiber cable models based upon literature
// Parameters:Fiber Diameter, number of segments, type, temperature
//	type: 1:Sundt Model 2:Tigerholm model 3:Rattay model
// All axon channels come from original models which were simplified to a 1D cable

public fiberD, nsegments, dx, Length,v
public section, sl, section_coord
objref section[1], section_coord, sl, conductances
create node[1]

proc global_parameters(){
	dx=segdensity //um//segment length which is for rattay which is most dense out of the fibers used
	nsegments = int(len/dx)
	Length = nsegments*dx
  }

proc build(){ 
  objref section[nsegments]
  section_coord = new Vector(nsegments,0)
  create node[nsegments]
  
  sl = new SectionList()
  for i=0,nsegments-1{
		node[i]{ 		
			ii = i
			section[i] = new SectionRef()
			sl.append(section[i])
			section_coord.x[ii] = (.5*dx + i*dx - Length/2) * 1e-6
			nseg = 1
			diam = fiberD
			L  = dx
			if ((i==0 || i==(nsegments-1)) && (type != 2)){
				insert pas
				g_pas = 0.0001
				if(type==1){
					e_pas=-60//Sundt equilibrium potential
				}else if(type == 2){
					e_pas = -55//tigerholm equilibrium potential
				}else if(type == 3){
					e_pas = -70//Rattay equilibrium potential
				}else if(type == 4){
					e_pas = -48//Sundt equilibrium potential
				}else{
					e_pas = -70
				}
				
				cm=1
				
				//no myelination
				insert extracellular
 				xg = 1e10 // short circuit, no myelin
				xc = 0    // short circuit, no myelin
				
				Ra = 1e10 // this forces current passive membrane and not into cable
			}else{
				if(type == 1){ // Sundt Fiber
					insert nahh
					gnabar_nahh = .04
					mshift_nahh = -6		        // NaV1.7/1.8 channelshift
					hshift_nahh = 6		        // NaV1.7/1.8 channelshift

					insert borgkdr			// insert delayed rectifier K channels
					gkdrbar_borgkdr = .04		// density of K channels
					ek = -90	  		// K equilibrium potential

					insert pas			// insert leak channels	
					g_pas = 1/10000			// set Rm = 10000 ohms-cm2
					Ra = 100			// intracellular resistance
					v=-60
					e_pas = v + (ina + ik)/g_pas	// calculate leak equilibrium potential
				}else if(type == 2){// Tigerholm Fiber
					conductances = new File()
					conductances.ropen("conductances.csv")
					insert ks
					gbar_ks 	= conductances.scanvar()
					insert kf
					gbar_kf 	= conductances.scanvar()
					insert h
					gbar_h 		= conductances.scanvar()
					insert nattxs
					gbar_nattxs 	= conductances.scanvar()
					insert nav1p8
					gbar_nav1p8 	= conductances.scanvar()
					insert nav1p9
					gbar_nav1p9 	= conductances.scanvar()
					insert nakpump
					smalla_nakpump 	= conductances.scanvar()
					insert kdrTiger
					gbar_kdrTiger 	= conductances.scanvar()
					insert kna
					gbar_kna 	= conductances.scanvar()
					insert naoi
					theta_naoi	= conductances.scanvar()	    
					insert koi
					theta_koi	= conductances.scanvar()	    
	
					insert leak
					insert extrapump	
	
					Ra 		= conductances.scanvar()
					cm 		= conductances.scanvar()
	
					conductances.close()
	
					celsiusT_ks 	= temp
					celsiusT_kf 	= temp
					celsiusT_h 	= temp
					celsiusT_nattxs = temp
					celsiusT_nav1p8 = temp
					celsiusT_nav1p9 = temp
					celsiusT_nakpump = temp
					celsiusT_kdrTiger = temp   
					v=-55
				}else if(type == 3){// Rattay and Aberham Fiber
					insert RattayAberham// Model adjusted for a resting potential of -70mV instead of 0 (subtract Vrest from each reversal potential)
					Ra=100
					cm=1
					v=-70
					ena=45
					ek=-82
				}else if(type == 4){//Schild Fiber
					R = 8314
					F = 96500
					insert leakSchild						// All mechanisms from Schild 1994 inserted into model
					insert kd
					insert ka
					insert can
					insert cat
					insert kds
					insert kca
					insert caextscale
					insert caintscale
					insert CaPump
					insert NaCaPump
					insert NaKpumpSchild
					if (insert97na) {
						insert naf97mean
						insert nas97mean
					} else {
						insert naf
						insert nas
					}
					// Ionic concentrations
					cao0_ca_ion = 2.0							// [mM] Initial Cao Concentration
					cai0_ca_ion = 0.000117						// [mM] Initial Cai Concentrations
					ko = 5.4									// [mM] External K Concentration
					ki = 145.0									// [mM] Internal K Concentration
					kstyle=ion_style("k_ion",1,2,0,0,0) 		// Allows ek to be calculated manually
					ek = ((R*(celsius+273.15))/F)*log(ko/ki)	// Manual Calculation of ek in order to use Schild F and R values
			
					nao = 154.0									// [mM] External Na Concentration
					nai = 8.9									// [mM] Internal Na Concentration
					nastyle=ion_style("na_ion",1,2,0,0,0) 		// Allows ena to be calculated manually
					ena = ((R*(celsius+273.15))/F)*log(nao/nai)	// Manual Calculation of ena in order to use Schild F and R values
					if (conductances97){				
						gbar_naf97mean=0.022434928 		// [S/cm^2] This block sets the conductance to the conductances in Schild 1997
						gbar_nas97mean=0.022434928
						gbar_kd=0.001956534
						gbar_ka=0.001304356
						gbar_kds=0.000782614
						gbar_kca=0.000913049
						gbar_can=0.000521743
						gbar_cat=0.00018261
						gbna_leakSchild = 1.8261E-05
						gbca_leakSchild	= 9.13049E-06
					}
					Ra=100
					cm = 1.326291192
					v=-48
					L_caintscale = L
					nseg_caintscale = nseg
					L_caextscale = L
					nseg_caextscale = nseg

				}
				insert extracellular
				xg = 1e10 // short circuit
				xc = 0    // short circuit
			}
		}
	}
	for i=0, nsegments-2 {
		connect node[i](1), node[i+1](0)	
	} 
}

proc init(){
  fiberD = 6 //default if no arg
  len = 21 //default if no arg
  type = 1
  temp = 37  //uF/cm2
  segdensity=50/6
  insert97na=0
  if (numarg()>0) {fiberD=$1}
  if (numarg()>1) {len=$2}
  if (numarg()>2) {type=$3}
  if (numarg()>3) {temp=$4}
  if (numarg()>4) {segdensity=$5}
  if (numarg()>5) {insert97na=$6}
  if (numarg()>6) {conductances97=$7}
  celsius = temp
  global_parameters() 
  build()
}

endtemplate cFiberBuilder