load_proc("nrnmainmenu")
load_file("nrngui.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(10000))+0.9)/2)*2+1}
}
	

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

proc parameters() {
	celsius = 33		
	Ri = 160	    			
	Cm = 1.0      		
	Rm = 15000      	
	v_init = -80
	Spinescale = 1.5					
}


// ----------------------------------------------------------------------------
// Channel densities
//-----------------------------------------------------------------------------

//CAVE: check for parameter-settings in loaded sessions
nagna=0.5

na_soma=220

na_121=250
na_122=900
na_123=350
na_124=225

na_ais1=0
na_ais2=500
na_ais3=700
na_ais4=500


na_node		=	7000
na_collat	=	500				
na_dend		=	20 						
na_myelin	=	40
max_distance_apic	=	300		// distance where gbar_na has reached na_dend
max_distance_basal	=	100		//corresponding distance for basal dendrites

vShift_na		=	10 // affecting activation and inactivation
vShift_inact_na		=	10 // affects only inactivation
q10_na			=	3
q10h_na			=	3

vShift_nax		=	10
vShift_inact_nax	=	10		
q10_nax			=	3.0
q10h_nax		=	3.0


length_constant_Kv_and_Kv1	=	80	//length constant of the exponential decay in um of Kv and Kv1
						//according to Keren et al., 2009, JPhysiol, 587:1413-37
Kv_soma	 	=	10
Kv_dend	 	=	3	
		

Kv1_dend	=	0.3		
Kv1_ais		=	1500 		
Kv1_soma	=	10
Kv1_collat	=	400

vShift_Kv1		=	10
vShift_inact_Kv1	=	-15

Kv7_soma	=	0.2	
Kv7_dend	=	0.2
Kv7_ais		=	0.300000
	
ca_reduce_fac	=	0.1
cagca=0.4

gca_dend	=	2.0*ca_reduce_fac
gca_soma	=	0.0005 
gca_ais1		=	0.0005
gca_ais2		=	3e-05 
gca_ais3		=	0

git2_soma	=	0.0003
git2_ais1	=	0.0001
git2_ais2	=	3e-05
git2_ais3	=	3e-06
git2_dend	=	2.0*ca_reduce_fac
git2_apic	=	6.0*ca_reduce_fac


gskca_soma	=	0.02
gbkca_soma	=	10
gskca_ais	=	0.02
gbkca_ais	=	1000  
gkca_dend	=	1.0*ca_reduce_fac

gCa_HVA_apic_hot_fac=	1	//i.e. no Ca hot spot	//ca
gCa_LVAst_apic_hot_fac=	1				//it2
gCa_hot_start	=	685
gCa_hot_end	=	885
vshiftm_cah=-4
vshifth_cah=-7.11157	
ct1_cah=0.1
ct2_cah=30

spinescale	=	1.5
sheaths		=	10	//number of myelin sheaths

cdelay1=1.8
cdelay2=315
cduration1=4
cduration2=7
camplitude1=-0.15
camplitude2=-0.3
cicin=0


// ----------------------------------------------------------------------------
// Initialisation of passive and active properties, spine scaling, Rm and Cm
//-----------------------------------------------------------------------------

proc init_channels() {

//all		
	forall {
		insert pas 	 
		Ra=Ri 
		e_pas=v_init
		g_pas=1/Rm
		cm=Cm
		insert nadifl
		D_nadifl= 0.4
	}

// dendrites
	for i=0,dendNum-1 {
		dend[i] {
		
			g_pas=1/(Rm/spinescale)
			cm=Cm*spinescale 
			
			insert na 	//gbar_na=na_dend*spinescale		//see below
			insert Kv 	gbar_Kv=Kv_dend*spinescale 		//see below
			insert Kv7 	gbar_Kv7=Kv7_dend
			insert kca 	gbar_kca = gkca_dend*spinescale 
			insert ca 	gbar_ca = gca_dend*spinescale
			insert it2 	gcabar_it2 = git2_dend*spinescale
			insert Kv1	gbar_Kv1 = Kv1_dend
			
		}			 
	}
	
	for i=0,apicNum-1 {
		apic[i] {
		
			g_pas=1/(Rm/spinescale)
			cm=Cm*spinescale
			insert na 	//gbar_na=na_dend*spinescale		//see below	
			insert Kv 	gbar_Kv=Kv_dend*spinescale 		//see below
			insert Kv7 	gbar_Kv7=Kv7_soma
			insert kca 	gbar_kca = gkca_dend*spinescale
			insert ca 	gbar_ca = gca_dend*spinescale
			insert it2 	gcabar_it2 = git2_apic*spinescale
			insert Kv1	gbar_Kv1 = Kv1_dend
		}			 
	}

//soma	
	soma  {nseg=9}
	soma  {
		insert na12 gbar_na12 = na_soma	

		insert Kv 	gbar_Kv = Kv_soma

		insert cah 	pbar_cah = gca_soma

		insert car	pbar_car = git2_soma

		insert bks	gkbar_bks = gbkca_soma

		insert Kv7 	gbar_Kv7 = Kv7_soma

		insert Kv1 	gbar_Kv1 = Kv1_soma

		insert cdp5

        insert cdp5r

        insert sk 	gbar_sk = gskca_soma

		//stimulation
		insert curclamp delay1_curclamp=cdelay1
		                delay2_curclamp=cdelay2
		                duration1_curclamp=cduration1
		                durataion2_curclamp=cduration2		               
		                amplitude1_curclamp=camplitude1		                
		                amplitude2_curclamp=camplitude2
		                icin_curclamp=cicin

	}

//collaterals note that axon[0] is AIS
	for i=1,axonNum-1 {	
		axon[i]  {
			g_pas=1/Rm
			insert nax gbar_nax = na_collat
			insert Kv1 gbar_Kv1 = Kv1_collat
			//insert cdp5
		}
	}	

//AIS
	axon[0]  {
		nseg=20

		insert na16 gbar_na16(0:nagna) = na_ais1:na_ais2
					gbar_na16(nagna:nagna+0.2) = na_ais2:na_ais3
					gbar_na16(nagna+0.2:1) = na_ais3:na_ais4
					
		insert na12 gbar_na12(0:0.1) = na_121:na_122
					gbar_na12(0.1:0.3) = na_122:na_122
					gbar_na12(0.3:0.5) = na_122:na_123
					gbar_na12(0.5:1) = na_123:na_124					

		insert Kv1 	axon.gbar_Kv1(0:1)= Kv1_soma:Kv1_ais

		insert Kv7 	axon.gbar_Kv7(0:0.4) = Kv7_soma:Kv7_soma
					axon.gbar_Kv7(0.4:1) = Kv7_soma:Kv7_ais
		//insert cad 
		insert cah 	pbar_cah (0:cagca)= gca_ais1:gca_ais2
					pbar_cah (cagca:1)= gca_ais2:gca_ais3

		insert car 	pbar_car (0:1)= git2_ais1:git2_ais2
					pbar_car (cagca:1)= git2_ais2:git2_ais3
		
		insert bks 	gkbar_bks(0:1) = gbkca_soma: gbkca_ais
	        	
		insert cdp5 TotalBuffer_cdp5 = 0.4
                    
        insert cdp5r 

        insert sk 	gbar_sk = gskca_ais

	}



//main axon: i.e. myelin and nodes
	for i=0,myNum-1 {
		my[i] {
			g_pas=1/(Rm * (sheaths + 1))
			cm=Cm/(sheaths + 1)			
			insert nax gbar_nax = na_myelin
			insert Kv7 gbar_Kv7 = Kv7_soma
			insert Kv1 gbar_Kv1 = 20				
		}	
	}	

	for i=0,myNum-1 {
		node[i]  {
			g_pas=1/Rm
			insert nax gbar_nax = na_node
			insert Kv1 gbar_Kv1 = Kv1_ais
			insert Kv7 gbar_Kv7 = Kv7_ais
		}
	}	


	forall ena=60
	forall ek=-88	

// ---------------------------------------------------------------------
// Calcium enhancement to reproduce frequency effect (Larkum et al,1999)
// ----------------------------------------------------------------------

	forall { 	
		vh1_it2    = 56
		vh2_it2    = 415
		ah_it2     = 30				
		v12m_it2   = 45
		v12h_it2   = 65  
		am_it2     = 3
		vshift_it2 = 10
		vm1_it2    = 50
		vm2_it2    = 125
	}

	forall if(ismembrane("ca_ion")) {
		eca = 140
		ion_style("ca_ion",0,1,0,0,0)
		vshift_ca = 8
	} 

	forall 	{
		caix_kca  = 0.7  	// Ca-sensitivity of Kca channel
		Ra_kca    = 1  	// Activation rate Kca 
		Rb_kca    = 2.5	// inactivation rate Kca, important ! 
        mi1_cah=1.2
        ti1_cah=1
	}

	//possibility to implement hot-zone
	access soma
	area(0.5)
	distance()
	//apical dendrites
	for i=0,apicNum-1 {
		apic[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				if(distance(pos) > gCa_hot_start && distance(pos) < gCa_hot_end) {
					gbar_ca(pos)	=	gca_dend*gCa_HVA_apic_hot_fac*spinescale
					gcabar_it2(pos)	=	git2_apic*gCa_LVAst_apic_hot_fac*spinescale
				}
			}
		}
	}

// --------------------------------------------------------------------------
// Dendritic exponential distribution of Ih (Kole et al., 2006)
// --------------------------------------------------------------------------
			 
	A = 0.0002 	//0.000429 before, and also the J Neurosci value I aim for a ~20 Mohm input resistance	
	tau = 0.003087
	spinescale = 1.5	

// --------------------------------------------------------------------------
// Dendritic linear distribution of na (Keren et al., 2009, JPhysiol, 587:1413-37)
// --------------------------------------------------------------------------

	access soma
	area(0.5)
	distance()
		
	//basal dendrites
	for i=0,dendNum-1 {
		dend[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				if (distance(pos) >= max_distance_basal)  {
					gbar_na(pos)=na_dend*spinescale	
				} else {
					gbar_na(pos)=(na_dend+(na_soma-na_dend)*(1-(distance(pos)/max_distance_basal)))*spinescale	
				}
			}
		}
	}

	//apical dendrites
	for i=0,apicNum-1 {
		apic[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				if (distance(pos) >= max_distance_apic)  {
					gbar_na(pos)=na_dend*spinescale	
				} else {
					gbar_na(pos)=(na_dend+(na_soma-na_dend)*(1-(distance(pos)/max_distance_apic)))*spinescale	
				}
			}
		}
	}




// --------------------------------------------------------------------------
// Dendritic linear distribution of Kv and Kv1 (Keren et al., 2009, JPhysiol, 587:1413-37)
// --------------------------------------------------------------------------

	access soma
	area(0.5)
	distance()
		
	//basal dendrites
	for i=0,dendNum-1 {
		dend[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				gbar_Kv(pos)=(Kv_dend+(Kv_soma-Kv_dend)*exp(-distance(pos)/length_constant_Kv_and_Kv1))*spinescale	
				gbar_Kv1(pos)=(Kv1_dend+(Kv1_soma-Kv1_dend)*exp(-distance(pos)/length_constant_Kv_and_Kv1))*spinescale	
			}
		}
	}

	//apical dendrites
	for i=0,apicNum-1 {
		apic[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				gbar_Kv(pos)=(Kv_dend+(Kv_soma-Kv_dend)*exp(-distance(pos)/length_constant_Kv_and_Kv1))*spinescale	
				gbar_Kv1(pos)=(Kv1_dend+(Kv1_soma-Kv1_dend)*exp(-distance(pos)/length_constant_Kv_and_Kv1))*spinescale	
			}
		}
	}
}


proc testAp() {
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "---------------------------------------------------------         A  P  I   C         ---------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"

	for i=0,apicNum-1 {
		apic[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				print " i ",i, " nseg ", j, " pos ",pos, " distance(pos) ",distance(pos)
				print " 		gbar_ih ", gbar_ih(pos), " gbar_na ", gbar_na(pos), " gbar_Kv ", gbar_Kv(pos), " gbar_Kv1 ", gbar_Kv1(pos)
			}
		}
	}
}


proc testBa() {
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "---------------------------------------------------------         B  A  S  A   L        -------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	print "-----------------------------------------------------------------------------------------------------------------------------------------------------------"
	for i=0,dendNum-1 {
		dend[i] {
			for j=1,nseg {
				pos=0.5/nseg+(j-1)/nseg
				print " i ",i, " nseg ", j, " pos ",pos, " distance(pos) ",distance(pos)
				print " 		gbar_ih ", gbar_ih(pos), " gbar_na ", gbar_na(pos), " gbar_Kv ", gbar_Kv(pos), " gbar_Kv1 ", gbar_Kv1(pos)
			}
		}
	}
}

access axon[0]
L=40
//undernai undersamples nai from 10kHz to 5kHz to match experiments
insert undernai
access soma