load_file("nrngui.hoc")

//The user can stimulate the terminals with capsaicin by adjusting Capsaicin_stimulation.hoc and uncommenting the last line of this file.

celsius = 25

//
// Nociceptive terminal model 
//
// Barkai et al. 2020 
// 

//////////////////// VARIABLES /////////////////////

tstop = 3000
CENTRAL_DIAM = 0.4

/////////////////   MORPHOLOGY ///////////////////////
Term_Branch = 27//Number of Terminal Branches
Term_Length = 1000 //Terminal Branches length
NavLessNum=14 
NavLessLength=25
SF=1 //Scaling factor for terminal 
				
				print NavLessLength
				CENTRAL_DIAM = 0.4
				create soma, central,terminal[Term_Branch], peri, tjcentral, tjperi, stem, cone, NavLess[NavLessNum]

				l_seg = 10 //segment size
				l_seg_skinterminal = 5 //segment size for terinals (for more accuracy)

				peri {nseg=100 L=5000 diam = .8}		// peripherial axon
				tjperi {nseg=100 L=100 diam = .8}		   	// tjunction axon
				tjcentral {nseg=100 L=100 diam = CENTRAL_DIAM}		// tjunction axon
				central {nseg=100 L=5000 diam = CENTRAL_DIAM}		// central axon
				stem {nseg=100 L=75 diam = 1.4}				// stem axon
				soma {nseg=1 L=25 diam =25}				// soma

				Rorg=1 //Mitochondira etc. resistance at terminals 

				//////// Connectivity //////////

				//terminal connect peri(0),1
				peri connect tjperi(0),1
				tjperi connect stem(0),1
				stem connect soma(0),1
				tjperi connect tjcentral(0),1
				tjcentral connect central(0),1

				//Terminal branch morphology and connectivity
				term_L=100
				term_nseg=term_L/l_seg   //Dividing the compartment to l size segment 
				term_diam=0.25

				cone connect peri(0),1
				connect terminal[0](1), cone(0)
				terminal[0]{nseg=term_nseg L=term_L diam=term_diam}
					
				// Terminal Morphology
					terminal[0] {nseg=term_nseg L=300*SF diam=term_diam}
					terminal[1] {nseg=term_nseg L=120*SF diam=term_diam}
					terminal[2] {nseg=term_nseg L=120*SF diam=term_diam}
					terminal[3] {nseg=term_nseg L=400*SF diam=term_diam}
					terminal[4] {nseg=term_nseg L=250*SF diam=term_diam}
					terminal[5] {nseg=term_nseg L=95*SF-NavLessLength diam=term_diam}
					terminal[6] {nseg=term_nseg L=175*SF-NavLessLength diam=term_diam} 
					terminal[7] {nseg=term_nseg L=40*SF diam=term_diam}
					terminal[8] {nseg=term_nseg L=40*SF diam=term_diam}
					terminal[9] {nseg=term_nseg L=40*SF diam=term_diam}
					terminal[10] {nseg=term_nseg L=160*SF diam=term_diam}
					terminal[11] {nseg=term_nseg L=100*SF-NavLessLength diam=term_diam}
					terminal[12] {nseg=term_nseg L=100*SF-NavLessLength diam=term_diam}
					terminal[13] {nseg=term_nseg L=100*SF-NavLessLength diam=term_diam}
					terminal[14] {nseg=term_nseg L=100*SF-NavLessLength diam=term_diam}
					terminal[15] {nseg=term_nseg L=350*SF diam=term_diam}
					terminal[16] {nseg=term_nseg L=175*SF-NavLessLength diam=term_diam}
					terminal[17] {nseg=term_nseg L=185*SF-NavLessLength diam=term_diam}
					terminal[18] {nseg=term_nseg L=85*SF diam=term_diam}
					terminal[19] {nseg=term_nseg L=200*SF-NavLessLength diam=term_diam}
					terminal[20] {nseg=term_nseg L=200*SF-NavLessLength diam=term_diam}
					terminal[21] {nseg=term_nseg L=155*SF-NavLessLength diam=term_diam}
					terminal[22] {nseg=term_nseg L=200*SF diam=term_diam}
					terminal[23] {nseg=term_nseg L=150*SF-NavLessLength diam=term_diam}
					terminal[24] {nseg=term_nseg L=80*SF diam=term_diam}
					terminal[25] {nseg=term_nseg L=75*SF-NavLessLength diam=term_diam}
					terminal[26] {nseg=term_nseg L=180*SF-NavLessLength diam=term_diam}
					
					
					connect terminal[1](1),terminal[0](0)	
					connect terminal[2](1),terminal[0](0) 
					connect terminal[3](1),terminal[1](0)	
					connect terminal[4](1),terminal[1](0) 
					connect terminal[5](1),terminal[2](0)
						connect NavLess[13](1), terminal[5](0)
					connect terminal[6](1),terminal[2](0)
						connect NavLess[12](1), terminal[6](0)
					connect terminal[7](1),terminal[3](0)	
					connect terminal[8](1),terminal[3](0)
					connect terminal[9](1),terminal[4](0)	
					connect terminal[10](1),terminal[4](0)	
					connect terminal[11](1),terminal[7](0)//skin
						connect NavLess[11](1), terminal[11](0)
					connect terminal[12](1),terminal[7](0)//skin
						connect NavLess[10](1), terminal[12](0)
					connect terminal[13](1),terminal[8](0)//skin	
						connect NavLess[9](1), terminal[13](0)
					connect terminal[14](1),terminal[8](0)//skin	
						connect NavLess[8](1), terminal[14](0)	
					connect terminal[15](1),terminal[9](0)	
					connect terminal[16](1),terminal[9](0)
						connect NavLess[7](1), terminal[16](0)
					connect terminal[17](1),terminal[10](0)
						connect NavLess[6](1), terminal[17](0)
					connect terminal[18](1),terminal[10](0)		
					connect terminal[19](1),terminal[15](0)//skin
						connect NavLess[5](1), terminal[19](0)	
					connect terminal[20](1),terminal[15](0)//skin
						connect NavLess[4](1), terminal[20](0)
					connect terminal[21](1),terminal[18](0)	
						connect NavLess[3](1), terminal[21](0)
					connect terminal[22](1),terminal[18](0)		
					connect terminal[23](1),terminal[22](0)//skin
						connect NavLess[2](1), terminal[23](0)	
					connect terminal[24](1),terminal[22](0)
					connect terminal[25](1),terminal[24](0)//skin
						connect NavLess[1](1), terminal[25](0)
					connect terminal[26](1),terminal[24](0)//skin
						connect NavLess[0](1), terminal[26](0)
					

					
				cone { //Cone between the terminal and the peripheral axon
					L=3 
					nseg=3 
					for i=0,nseg{ //coning the cone with respect to the terminal and peri diameters
							
							diam((1/nseg)*i/2)=(1/nseg)*i*(peri.diam - terminal[0].diam) + terminal[0].diam    
					}
					diam(1)=peri.diam
				}

				objref DRG   					// create new object called DRG
				DRG = new SectionList()				// Define DRG as a list of sections

				tjperi DRG.append()				// DRG sections: tjunction, stem, and soma   
				tjcentral DRG.append()  
				stem DRG.append()				
				soma DRG.append()

				objref TermSec						// create new object called TermSec for terminals
				TermSec = new SectionList()	  					
					cone TermSec.append()
					for i=0, Term_Branch-1 { 	
							terminal[i] TermSec.append()
					}
					
				objref NavLessSec
				NavLessSec = new SectionList()	  					
					for i=0, NavLessNum-1 { 	
							NavLess[i] NavLessSec.append()
					}
					
				objref TermSkin	                    // Skin inervating Terinals (for Goldstein et al.)
				TermSkin = new SectionList()	
					terminal[26] TermSkin.append()
					terminal[25] TermSkin.append()
					terminal[23] TermSkin.append()
					terminal[21] TermSkin.append()
					terminal[20] TermSkin.append()
					terminal[19] TermSkin.append()
					terminal[17] TermSkin.append()
					terminal[16] TermSkin.append()
					terminal[14] TermSkin.append()
					terminal[13] TermSkin.append()
					terminal[12] TermSkin.append()
					terminal[11] TermSkin.append()
					terminal[5] TermSkin.append()
					terminal[6] TermSkin.append()
					



				/////////////////// INITIALIZATION /////////////////
				 

				///// GENERAL INITIALIZATION /////////

				forall{				// for all compartments
						insert na_ion  //A bug fix for point_process 
						insert k_ion   //A bug fix for point_process
					   
					   //Na Channels 
					   insert nattxs
							gbar_nattxs=0.0017
					   insert nav1p9
							gbar_nav1p9=0.00064
					   insert nav1p8
							gbar_nav1p8=0.02
						insert nap
							gbar_nap=5e-5
					   
					   
					   //K Channels
					   insert kdr
							gbar_kdr=0.00083
					   insert ka
							gkbar_ka=0.0015
					   insert km
							gbar_km=0.00034
					   //insert kf
					   
						//Ca channels
						insert calL_Shah
							gcalbar_calL=0.003
						insert calT_Shah 
							gcatbar_calT_Shah=0.001
						//H channels
						insert hd	
							ghdbar_hd=0.00033

						
					insert pas			// insert leak channels	
					g_pas =1/10000	 //set Rm = 10000 ohms-cm2
					v = -60				// set Vrest
					e_pas = -60
					Ra = 150			// intracellular resistance
				}


				forsec NavLessSec{
						Ra=15*peri.Ra
						g_pas=g_pas/4          //Accoring to  Dmytro V. Vasylyev and Stephen G. Waxman, Journal of Neurophysiology 2012	
						
						nseg=L/l_seg
						L=NavLessLength //-ShortenBy
						
						
						
						
						diam=term_diam
						
						distance() 
						
						insert transducer_pas					
						//Distrubuting No NaV channels in the distal end of the terminal 
											
						gbar_nattxs=0
						gbar_nav1p9=0
						gbar_nav1p8=0
						gbar_nap=0
										
						g_transMAX=0.0025
						Tau_Puff=7.04             //(from Golstein et al. 2017)
						for (x){ 
								xdist = distance(x)
								g_transducer_pas(x)=g_transMAX*exp(-(x*L)/Tau_Puff) //look at Capsaicin_Puff_Exponent_Check.m if you want to verify this
								}
				}	
				forsec TermSec{
					g_pas=g_pas/4          //Accoring to  Dmytro V. Vasylyev and Stephen G. Waxman, Journal of Neurophysiology 2012	
				}

				forsec TermSkin{
					
					//L=L+ShortenBy
					//nseg=L/l_seg
					
					insert transducer_pas
					Ra=15*peri.Ra	
					g_transMAX=0.0025	
					Tau_Puff=7.04 	
					g_MIN=g_transMAX*exp(-(NavLessLength)/Tau_Puff)
								//(from Golstein et al. 2017)
					distance()
					
					nseg=L/l_seg_skinterminal //make the terminal of stimulation with a higher (nseg) resolution	
					
					for (x){ 
							xdist = distance(x)				
							g_transducer_pas(x)=g_MIN*exp(-(x*L)/Tau_Puff) //look at Capsaicin_Puff_Exponent_Check.m if you want to verify this	
							}
							
				}


				proc init(){				// INITIALIZATION FUNCTION
						dt=0.5 //0.1	//dt=0.7 

				  forall {
					   v=-60 // VREST FOR ALL COMPARTMENTS (Change to -54mV to get rid of the spontaneous first spike)
					   ena=60
					   ek=-85
					   ehd_hd=-20
					   
					 
					   finitialize(v)			// reset all state variables
					   fcurrent()     			// calculate all currents
					   }
				}			/// end of initialization
		

			//The user can also choose which terminals to stimulate by changing Capsaicin_stimulation.hoc file and uncommenting the next line
			//load_file("Capsaicin_stimulation.hoc") // puts the electrodes in specific or all terminals