/*--------------------------------------------------------------------
2/02
Cameron C. McIntyre
SIMULATION OF PNS MYELINATED AXON

This model is described in detail in:

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.

This model can not be used with NEURON v5.1 as errors in the
extracellular mechanism of v5.1 exist related to xc. The original
stimulations were run on v4.3.1. NEURON v5.2 has corrected the 
limitations in v5.1 and can be used to run this model.
----------------------------------------------------------------------*/

load_proc("nrnmainmenu")

proc model_globels() {			
	celsius=37			
	v_init=-80 //mV//  		
	dt=0.005 //ms//         	
	tstop=10
//Intracellular stimuluation parameters//	
	istim=2				
	delay=1				
	pw=0.1
//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() {
	if (fiberD==5.7) {g=0.605 axonD=3.4 nodeD=1.9 paraD1=1.9 paraD2=3.4 deltax=500 paralength2=35 nl=80}
	if (fiberD==7.3) {g=0.630 axonD=4.6 nodeD=2.4 paraD1=2.4 paraD2=4.6 deltax=750 paralength2=38 nl=100}
	if (fiberD==8.7) {g=0.661 axonD=5.8 nodeD=2.8 paraD1=2.8 paraD2=5.8 deltax=1000 paralength2=40 nl=110}
	if (fiberD==10.0) {g=0.690 axonD=6.9 nodeD=3.3 paraD1=3.3 paraD2=6.9 deltax=1150 paralength2=46 nl=120}
	if (fiberD==11.5) {g=0.700 axonD=8.1 nodeD=3.7 paraD1=3.7 paraD2=8.1 deltax=1250 paralength2=50 nl=130}
	if (fiberD==12.8) {g=0.719 axonD=9.2 nodeD=4.2 paraD1=4.2 paraD2=9.2 deltax=1350 paralength2=54 nl=135}
	if (fiberD==14.0) {g=0.739 axonD=10.4 nodeD=4.7 paraD1=4.7 paraD2=10.4 deltax=1400 paralength2=56 nl=140}
	if (fiberD==15.0) {g=0.767 axonD=11.5 nodeD=5.0 paraD1=5.0 paraD2=11.5 deltax=1450 paralength2=58 nl=145}
	if (fiberD==16.0) {g=0.791 axonD=12.7 nodeD=5.5 paraD1=5.5 paraD2=12.7 deltax=1500 paralength2=60 nl=150}
	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()

objectvar stim

create node[axonnodes], MYSA[paranodes1], FLUT[paranodes2], STIN[axoninter]
access node[0]	//APD

proc initialize(){
	for i=0,axonnodes-1 {
		node[i]{					
			nseg=1
			diam=nodeD
			L=nodelength
			Ra=rhoa/10000
			cm=2
			insert axnode			
			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=-80
			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=-80
			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=-80
			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()

//intracellular stimulus//

node[10] stim = new IClamp(0.5)

proc setstimparams() {
  stim.del = delay
  stim.dur = pw
  stim.amp = istim
}

setstimparams()

xpanel("Stimulus parameters")
	xvalue("Stimulus Amplitude (nA)", "istim", 1, "setstimparams()", 1)
	xvalue("Pulse Duration (ms)", "pw", 1, "setstimparams()", 1)
	xvalue("Onset Delay (ms)", "delay", 1, "setstimparams()", 1)
xpanel(100,100)