/*--------------------------------------------------------------------
4/2013
Pietro Balbi
SIMULATION OF PNS MYELINATED AXON

Taken from:
McIntyre CC, Richardson AG, and Grill WM. Modeling the excitability of
mammalian nerve fibers: influence of afterpotentials on the recovery
cycle. Journal of Neurophysiology 87:995-1006, 2002.

----------------------------------------------------------------------*/

proc model_globels() {			
//topological parameters//		
	axonnodes=21  			
	paranodes1=40
	paranodes2=40	
	axoninter=120			
	axontotal=221			 
//morphological parameters//	
	fiberD=10.0	//choose from 5.7, 7.3, 8.7, 10.0, 11.5, 12.8, 14.0, 15.0, 16.0
	paralength1=3  
	nodelength=1.0
	space_p1=0.002  
	space_p2=0.004
	space_i=0.004
//electrical parameters//		
	rhoa=0.7e6 //Ohm-um//
	mycm=0.1 //uF/cm2/lamella membrane//
	mygm=0.001 //S/cm2/lamella membrane//
	}
model_globels ()

proc dependent_var() {  // values for fiberD=10
	g=0.690
	axonD=6.9
	nodeD=3.3
	paraD1=3.3
	paraD2=6.9
	deltax=1150
	paralength2=46
	nl=120
	Rpn0=(rhoa*.01)/(PI*((((nodeD/2)+space_p1)^2)-((nodeD/2)^2)))
	Rpn1=(rhoa*.01)/(PI*((((paraD1/2)+space_p1)^2)-((paraD1/2)^2)))
	Rpn2=(rhoa*.01)/(PI*((((paraD2/2)+space_p2)^2)-((paraD2/2)^2)))
	Rpx=(rhoa*.01)/(PI*((((axonD/2)+space_i)^2)-((axonD/2)^2)))
	interlength=(deltax-nodelength-(2*paralength1)-(2*paralength2))/6
	}
dependent_var()

create node[axonnodes], MYSA[paranodes1], FLUT[paranodes2], STIN[axoninter]

proc initialize(){
	for i=0,axonnodes-1 {
		node[i]{					
			nseg=1
			diam=nodeD
			L=nodelength
			Ra=rhoa/10000
			cm=2
			insert pas
				g_pas=0.001
				e_pas=-72
			insert na3rp
				gbar_na3rp=0.88
				sh_na3rp=5
				ar_na3rp=0.4
			insert naps
				gbar_naps=0.0044
				sh_naps=15
				ar_naps=0.4
			insert kdrRL
				gMax_kdrRL=0.4
			insert extracellular xraxial=Rpn0 xg=1e10 xc=0
			}
	}
	for i=0, paranodes1-1 {
		MYSA[i]{
			nseg=1
			diam=fiberD
			L=paralength1
			Ra=rhoa*(1/(paraD1/fiberD)^2)/10000
			cm=2*paraD1/fiberD
			insert pas
				g_pas=0.001*paraD1/fiberD		
				e_pas=v_init
			insert extracellular xraxial=Rpn1 xg=mygm/(nl*2) xc=mycm/(nl*2)
			}
	}
	for i=0, paranodes2-1 {
		FLUT[i]{
			nseg=1
			diam=fiberD
			L=paralength2
			Ra=rhoa*(1/(paraD2/fiberD)^2)/10000
			cm=2*paraD2/fiberD
			insert pas
				g_pas=0.0001*paraD2/fiberD		
				e_pas=v_init
			insert extracellular xraxial=Rpn2 xg=mygm/(nl*2) xc=mycm/(nl*2)
			}
	}
	for i=0, axoninter-1 {
		STIN[i]{
			nseg=1
			diam=fiberD
			L=interlength
			Ra=rhoa*(1/(axonD/fiberD)^2)/10000
			cm=2*axonD/fiberD
			insert pas
				g_pas=0.0001*axonD/fiberD
				e_pas=v_init
			insert extracellular xraxial=Rpx xg=mygm/(nl*2) xc=mycm/(nl*2)
			}
	}
	for i=0, axonnodes-2 {
		connect MYSA[2*i](0), node[i](1)
		connect FLUT[2*i](0), MYSA[2*i](1)
		connect STIN[6*i](0), FLUT[2*i](1)
		connect STIN[6*i+1](0), STIN[6*i](1)
		connect STIN[6*i+2](0), STIN[6*i+1](1)
		connect STIN[6*i+3](0), STIN[6*i+2](1)
		connect STIN[6*i+4](0), STIN[6*i+3](1)	
		connect STIN[6*i+5](0), STIN[6*i+4](1)	
		connect FLUT[2*i+1](0), STIN[6*i+5](1)
		connect MYSA[2*i+1](0), FLUT[2*i+1](1)
		connect node[i+1](0), MYSA[2*i+1](1)	
		}
	
	finitialize(v_init)
	fcurrent()
}
initialize()