xopen("A3+AIS_original.hoc")

access soma
area(0.5)
distance()

proc geom_nseg() {

soma area(0.5) // make sure diam reflects 3d points
forall {nseg = int((L/(0.1*lambda_f(2000))+0.9)/2)*2+1}

}
	



//--------------------------------------------------------------
//  Initialise user-defined membrane parameters 
// --------------------------------------------------------------

	proc parameters() {

	celsius=37		// nominal temperature of simulation
	Ri=100    		// intracellular resistivity ohm-cm
	Cm=0.9      		// specific membrane capacitance in uF/cm^2
	Rm=15000      		// specific membrane resistivity in ohm/cm^2 
	v_init=-88     	// to obtain a membrane potential of -75 at the soma 
	spinescale=2		// factor used to membrane area to account for spines
					
}



// --------------------------------------------------------------------
//   Axon and collaterals. Properties based on Sloper and Powell, 1979  
// --------------------------------------------------------------------	


create myelin[45], node[42]


	for i=0,44	{
    			myelin[i] {
      			nseg = 20	//for space plot
      			L = 60	//spacing 100 times the diameter  
      			diam = 1.6 
				}
    	}

	for i=0,41	{

			node[i] {
      			nseg = 1
      			L = 1           
      			diam = 1.1
				}
	}



//------------------------------------
//	Axon topology 
//-------------------------------------

	axon connect myelin[0](1), 1
	myelin[0] connect node[0](1), 0

	for i=0,8  {node[i] connect myelin[i+1](1), 0 
      		 myelin[i+1] connect node[i+1](1), 0
}


	node[1] connect myelin[10](1), 0
	myelin[10] connect node[10](1), 0

	for i=10,13  {node[i] connect myelin[i+1](1), 0 
      		myelin[i+1] connect node[i+1](1), 0
}


	node[0] connect myelin[15](1), 0
	myelin[15] connect node[15](1), 0
	node [15] connect myelin[44](1), 0

	node[2] connect myelin[42](1), 0
	myelin [42] connect myelin[43](1), 0

	node[12] connect myelin[16](1), 0
	myelin[16] connect node[16](1), 0

	for i=16,29  {node[i] connect myelin[i+1](1), 0 
      	 	myelin[i+1] connect node[i+1](1), 0
}

	node [20] connect myelin[31](1), 0
	myelin [31] connect node[31](1), 0

	for i=31,34  {node[i] connect myelin[i+1](1), 0 
      		myelin[i+1] connect node[i+1](1), 0
}

	node[33] connect myelin[36](1), 0
	myelin [36] connect node[36](1), 0

	node[26] connect myelin[37](1), 0
	myelin [37] connect node[37](1), 0

	for i=37,39  {node[i] connect myelin[i+1](1), 0 
      		 myelin[i+1] connect node[i+1](1), 0
}

	node[16] connect myelin[41](1), 0
	myelin [41] connect node[41](1), 0






// ----------------------------------------------------------------------------
// initialize passive and active properties and add spines by scaling Rm and Cm
//-----------------------------------------------------------------------------


na_segment = 3000
na_soma = 60
na_node= 2500 


proc init_channels() {

	na_myelin=80			// in dendrites pS/um2
						
	Km_soma=5				// units in pS/um2
	Km_axon=50			

	Kv_soma=20			// (HVA Kv current, units in pS/um2)
	Kv_axon=2000				

	Kv1_soma=0.01			// 0.01 mOhm/cm2 = 100 pS/um2 (LVA-Kv current)
	Kv1_axon=0.20 			
	
	spinescale=2
	
	vshift_na=10			//provides AP threshold of ~ -60 mV at soma 
	vshift_nax=10	
	


// dendrites

		forall {insert pas g_pas=1/(Rm/spinescale) 
				cm=Cm*spinescale 
				Ra=Ri 					
				e_pas=v_init
					}

		forall {insert Kv1 gbar_Kv1=Kv1_soma*spinescale
				insert Km gbar_Km=Km_soma*spinescale
				insert Kv gbar_Kv=Kv_soma*spinescale
					}		

		for i=0,42 {
				dend[i] {insert na gbar_na=na_soma*spinescale}			 
					}
		
		for i=0,86 {
				apic[i] {insert na gbar_na=na_soma*spinescale}			 
					}



// soma-hillock area

		soma  {nseg=8}
		soma  {g_pas=1/Rm cm=Cm}
		soma  {insert na gbar_na = na_soma} 
		soma  {insert Kv1 gbar_Kv1 = Kv1_soma} 
		soma  {insert Kv gbar_Kv = Kv_soma}	
		soma  {insert Km gbar_Km = Km_soma} 		
		
	
				

// Subcellular properties of initial segment
		
		axon  {nseg=20}
		axon	 {Rm=15000}		
		axon  {g_pas=1/Rm cm=Cm}
		axon	 {Ri=100}
		axon  {e_pas=-75}
		axon  {insert nax axon.gbar_nax(0:0.725)= na_segment:na_segment
				   axon.gbar_nax(0.775:1)= na_segment:na_soma}
		
		axon {insert Kv1 axon.gbar_Kv1(0:0.725) = Kv1_soma:Kv1_axon
					axon.gbar_Kv1(0.775:1) = Kv1_axon:Kv1_axon}
		
		axon {insert Kv axon.gbar_Kv(0:0.725) = Kv_soma:Kv_axon
					axon.gbar_Kv(0.775:1) = Kv_axon:Kv_axon}	

		axon {insert Km axon.gbar_Km(0:0.725) = Km_soma:Km_axon
					axon.gbar_Km(0.775:1) = Km_axon:Km_axon}
		
			

// Myelinated axon properties


	for i=0,44 {
		myelin[i] {g_pas=1/Rm  cm=0.02} 
		myelin[i] {e_pas=-75}			
		myelin[i] {insert Km gbar_Km = Km_soma}
		myelin[i]	{insert Kv gbar_Kv = Kv_soma}
		myelin[i] {insert na gbar_na = na_myelin}   
	}	
	

	for i=0,41 {
		node[i]  	{g_pas=400/Rm cm=Cm}
		node[i]	{e_pas=-75}	
		node[i]  	{insert nax gbar_nax = na_node} 
		node[i]  	{insert Km gbar_Km = Km_axon}
		node[i]  	{insert Kv1 gbar_Kv1 = Kv1_axon} 
		node[i]	{insert Kv gbar_Kv = Kv_axon} 		
	}

	forall ena=55
	forall ek=-85
	
}
	



// --------------------------------------------------------------------------
// Dendritic exponential distribution of Ih (Kole et al., 2006)
// --------------------------------------------------------------------------
		 

	hfactor = 0.000429 		
	hinvtau = 0.003087
	spinescale = 2	

	forall {insert ih
			gh_ih = -0.0002
			ehd_ih=-45		 
	
		} 
	
	access soma
	area(0.5)
	distance()


//exponential in basal dendrites
	for i=0,42 {
		access dend[i]
		ghdbar_ih=gh_ih+(hfactor*(exp(hinvtau*distance(.5))))
		ghdbar_ih=ghdbar_ih*spinescale	
	}


//exponential in apical dendrites

	for i=0,86 {
		access apic[i]
		ghdbar_ih=gh_ih+(hfactor*(exp(hinvtau*distance(.5))))
		ghdbar_ih=ghdbar_ih*spinescale	
	}

	soma {ghdbar_ih=ghdbar_ih/spinescale}	



// remove ih from the reconstructed axon 
	
	axon {ghdbar_ih=0}
		

	for i=0,44 {
		myelin[i] {ghdbar_ih=0} 
		}
	for i=0,41 {
		node[i]  	{ghdbar_ih=0}

		}