Chirp stimulus responses in a morphologically realistic model (Narayanan and Johnston, 2007)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:147578
...we built a multicompartmental model with a morphologically realistic three-dimensional reconstruction of a CA1 pyramidal neuron. The only active conductance we added to the model was the h conductance. ... We conclude that experimentally observed gradient in density of h channels could theoretically account for experimentally observed gradient in resonance properties (Narayanan and Johnston, 2007).
Reference:
1 . Narayanan R, Johnston D (2007) Long-term potentiation in rat hippocampal neurons is accompanied by spatially widespread changes in intrinsic oscillatory dynamics and excitability. Neuron 56:1061-75 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Dendrite;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampus CA1 pyramidal GLU cell;
Channel(s): I h;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Ion Channel Kinetics; Active Dendrites; Detailed Neuronal Models;
Implementer(s): Narayanan, Rishikesh [rishi at iisc.ac.in];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal GLU cell; I h;
/
MultiChirp
Output
README.txt
h.mod *
ChirpI.txt *
Fig8C.hoc
mosinit.hoc
n123.hoc *
ObliquePath.hoc *
oblique-paths.hoc *
Trunk_RN.txt
Trunk_Z.txt
                            
/**************************************************************
Written by Rishikesh Narayanan.

Papers:

Rishikesh Narayanan and Daniel Johnston, “Long-term potentiation
in rat hippocampal neurons is accompanied by spatially widespread
changes in intrinsic oscillatory dynamics and excitability,” Neuron,
56(6), pp. 1061–1075, December 2007.

Rishikesh Narayanan and Daniel Johnston, “The h channel mediates
location dependence and plasticity of intrinsic phase response in
rat hippocampal neurons,” The Journal of Neuroscience, 28(22), pp.
5846–5860, May 2008.
**************************************************************/


load_file("stdlib.hoc")
load_file("ObliquePath.hoc")

//**************************************************************//
// Current clamp Parameters
//**************************************************************//

tstop = 1500
steps_per_ms = 40
dt = 0.025

stamp=-0.1		// Stimulus Amplitude and Delay
stdur=1000
stdel=20

//**************************************************************//
// Chirp Parameters
//**************************************************************//

objref fchirp, ChirpStim
tstop = 25050
steps_per_ms = 40
dt = 0.025
ChirpAmp=0.05

// --------------------------------------------------------------
// Resistances and Capacitances
// --------------------------------------------------------------

rasoma  = 50			
raend   = 35
raaxon = 50

rm     = 200000
rmsoma = rm
rmdend = rm
rmend=12000
rmaxon = rm

c_m       = 1
cmsoma    = c_m
cmdend    = c_m
cmaxon    = c_m

v_init    = -65
celsius   = 34

// --------------------------------------------------------------
// Reversal potentials
// --------------------------------------------------------------

Eh=-30

//--------------------------------------------------------------
// Active conductance densities and other related parameters
//--------------------------------------------------------------

gh=0.0006		// Base value of gh

gh=3.4e-5

//--------------------------------------------------------------

create axon[2]
objref sh, st, apc, vec, g
strdef str1, str2, str3
objref SD, AXON, SA, Basal, Trunk
objref s, ampasyn[1], nmdasyn[1]
objref bpropampasyn, bpropnmdasyn
objref LM, RAD, RADt, LUC, PSA, PSB, ORI
objref secref, f1
objref netlist
objref apamp[1], ampvec[1]
objref sapamp, sampvec
objref g1, nc
objref r, st1
objref Trial
objref local_trunk
objref pl[150],opl[150]

create soma[1], apical[1], basal[1]

objref v3


/**********************************************************************/
// Passive Conductances 

proc setpassive() {
	forall {
	  	insert pas
	  	e_pas = v_init
		Ra=rasoma
	}
	forsec SD {		// For somato-dendritic compartments
		cm=cmdend
		g_pas=1/rmdend
	}
	forsec "soma" {
	  	cm = cmsoma 
		g_pas=1/rmsoma
	}
	forsec Trunk {
		for (x) {
			rdist=raddist(x)
			rm = rmsoma + (rmend - rmsoma)/(1.0 + exp((200-rdist)/50))
			g_pas(x)=1/rm
			Ra = rasoma + (raend - rasoma)/(1.0 + exp((210-rdist)/50))
		}
	}	

	for i=0,plcount {
    	seccount=0
    	forsec pl[i] {
       		if(!seccount){
           		trunk_pas=g_pas(1)
				seccount=seccount+1
       		} else {
           		g_pas=trunk_pas
				seccount=seccount+1
       		}
       		//print secname() 
    	}
	}       

}

/**********************************************************************/
// Active Conductances 
// Ih has somatic conductance, as in Poirazi et al., 2003.

proc setactive () {local xdist, ndist

	forsec Trunk {	// Trunk
		insert hd  

		for (x) {
			xdist=raddist(x)
			ghdbar_hd(x) = gh*(1+100/(1+exp(-(xdist-380)/34.0)))

			if (xdist > 100){
				if (xdist>300) { 
					ndist=300
				} else { // 100 <= xdist <= 300
					ndist=xdist
				}
				vhalfl_hd(x)=-82-8*(ndist-100)/200
			} else {	// xdist < 100
				vhalfl_hd(x)=-82
			}	

		}
	}

	for i=0,plcount { // Aplical obliques
    	seccount=0
    	forsec pl[i] {
        	if(!seccount){	// The first section is the trunk
            	trunk_h=ghdbar_hd(1)
				trunk_vhalf=vhalfl_hd(1)
				seccount=seccount+1
        	} else {
				insert hd
            	ghdbar_hd=trunk_h
				vhalfl_hd=trunk_vhalf
				seccount=seccount+1
        	}
        	//print secname() 
    	}
	}       

	forsec "soma" {
		insert hd  ghdbar_hd=gh vhalfl_hd=-82
	}	

	forsec Basal {
		insert hd  ghdbar_hd=gh vhalfl_hd=-82
	}	

	forall if (ismembrane("hd") ) ehd_hd=Eh
	
}

/**********************************************************************/

proc init_cell() {
	setpassive()
	addaxon()
	setactive()
	access soma[2]	// Reinitializing distance origin
	distance()
	print "Number of compartments: ", totcomp


	finitialize(v_init)
	fcurrent()
        forall {
            for (x) {
        		if (ismembrane("hd")) {
					e_pas(x)=v(x)+i_hd(x)/g_pas(x)
                } else {
        			e_pas(x)=v(x)
                }
            }
        }
}

/**********************************************************************/

proc current_clamp(){	
	access apical[65] 
	st1=new IClamp(.5)
	st1.dur = stdur
	st1.del = stdel
	st1.amp = stamp
	access soma[0]
}

/**********************************************************************/

proc load_3dcell() {

// $s1 filename

	forall delete_section()
	xopen($s1)

	SD = new SectionList()
	SA = new SectionList()
	Trunk = new SectionList()
	Trial = new SectionList()
	Basal = new SectionList()

	forsec "soma" {
		SD.append()
		SA.append()
	}

	forsec "basal" { 
		SD.append()
		Basal.append()
	}

	forsec "apical"{
		SD.append()
		SA.append()
	}

	// Trunk. 

	soma[0]	   Trunk.append()
	apical[0]  Trunk.append() 
	apical[4]  Trunk.append() 
	apical[6]  Trunk.append() 
	apical[14] Trunk.append() 
	apical[15] Trunk.append() 
	apical[16] Trunk.append() 
	apical[22] Trunk.append() 
	apical[23] Trunk.append() 
	apical[25] Trunk.append() 
	apical[26] Trunk.append() 
	apical[27] Trunk.append() 
	apical[41] Trunk.append() 
	apical[42] Trunk.append() 
	apical[46] Trunk.append() 
	apical[48] Trunk.append() 
	apical[56] Trunk.append() 
	apical[58] Trunk.append() 
	apical[60] Trunk.append() 
	apical[62] Trunk.append() 
	apical[64] Trunk.append() 
	apical[65] Trunk.append() 
	apical[69] Trunk.append() 
	apical[71] Trunk.append()
	apical[81] Trunk.append() 
	apical[83] Trunk.append() 
	apical[95] Trunk.append()
	apical[103] Trunk.append()
	apical[104] Trunk.append()

	// Define obliques

	load_file("oblique-paths.hoc")

	// Compartmentalize 

	setpassive() // Before compartmentalization, set passive

	// The lambda constraint
	totcomp=0
	forall { 
		soma[1] area(0.5)
		nseg = int((L/(0.1*lambda_f(100))+.9)/2)*2 + 1  
		totcomp=totcomp+nseg
	}
	access soma[2]
	distance()

	init_cell()

// If you want to know the Radial distances of all trunk compartments.

	/*forsec Trunk {
		for(x) {
			print secname(), " ", x, " ", raddist(x)
		}
	}*/	
}

/**********************************************************************/
// For cell number n123 on the DSArchive, converted with CVAPP to give
// HOC file, the following definition holds. This is the same as Poirazi et
// al. have used in Neuron, 2003. The argument is that the subtree seems
// so long to be a dendrite, and the cell does not have a specific axon.
// There is a catch, though, if the morphology is closely scanned, then 
// basal dendrites would branch from these axonal segments - which
// may be fine given the amount of ambiguity one has while tracing! 

proc addaxon() {

	AXON = new SectionList()

	for i = 30,34 basal[i] {
		AXON.append()
		Basal.remove()
	}

	for i = 18,22 basal[i] {
		AXON.append()
		Basal.remove()
	}

	forsec AXON {
		e_pas=v_init 
		g_pas = 1/rmaxon 
		Ra=raaxon 
		cm=cmaxon
	}
}

/********************************************************************/

proc initsub() {
	load_3dcell("n123.hoc")
	finitialize(v_init)
	fcurrent()
}

/********************************************************************/

proc update_Gh() { local xdist
	forsec "apical" {
	   if (ismembrane("hd") ) {
			for (x) {
				xdist=raddist(x)
				ghdbar_hd(x)=gh*(1+2*xdist/100.0)
			}
	    }
	}

	forsec "soma" {
		ghdbar_hd=gh 
	}	

	forsec Basal {
		ghdbar_hd=gh
	}	

	update_init()
}

/********************************************************************/
proc update_init(){
	finitialize(v_init)
	fcurrent()
        forall {
            for (x) {
        		if (ismembrane("hd")) {
					e_pas(x)=v(x)+i_hd(x)/g_pas(x)
                } else {
        			e_pas(x)=v(x)
                }
            }
        }
}
/********************************************************************/
// Procedure to calculate radial distances. The following three coordinates
// are the centroid of the somatic compartments. Radial distances are
// computed by translating each x to corresponding pt3d()s and using this
// centroid for computing the radial distance.

somax=2.497
somay=-13.006
somaz=11.12
double distances[200]

func raddist() {
	distn0=distance(0)
	distances[0]=0
	sum=0

	for i=1,n3d()-1 {
		xx=(x3d(i)-x3d(i-1))*(x3d(i)-x3d(i-1))
		yy=(y3d(i)-y3d(i-1))*(y3d(i)-y3d(i-1))
		zz=(z3d(i)-z3d(i-1))*(z3d(i)-z3d(i-1))
		sum=sum+sqrt(xx+yy+zz)
		distances[i]=sum
	}

	xval=$1

// Amoung the various pt3d's find which one matches the distance of
// current x closely

	distn=distance(xval)
	match=distn-distn0
	matchptdist=100000
	for i=0,n3d()-1 {		
		matptdist=(match-distances[i])*(match-distances[i])
		if(matchptdist>matptdist){
			matchptdist=matptdist
			matchi=i
		}
	}

// Find the distance of the closely matched point to the somatic
// centroid and use that as the distance for this BPAP measurement			

	xx=(x3d(matchi)-somax)*(x3d(matchi)-somax)
	yy=(y3d(matchi)-somay)*(y3d(matchi)-somay)
	zz=(z3d(matchi)-somaz)*(z3d(matchi)-somaz)
	return sqrt(xx+yy+zz)
}

/*****************************************************************/

proc GenerateChirp() {
	ChirpStim=new Vector(25050/dt)
	ChirpStim.fill(0)
	for i=50/dt, ChirpStim.size()-1 {
		xval=(i-50/dt)*dt/1000
		ChirpStim.x[i]=ChirpAmp*sin(2*3.141592654*25*(i-50/dt)/(2*(25000/dt)) *xval)
	}
	fchirp=new File()
	fchirp.wopen("ChirpI.txt")
	ChirpStim.printf(fchirp,"%f\n")
}


/*****************************************************************/
// Compute the response of the trunk compartments to local Chirp stimuli
// Output saved in files in a directory called Output. Create a directory
// called Output before you run this.

proc Chirp_Trunk() {
	print "\n\nComputing Trunk Compartment Responses to Local Chirp Stimuli.."
	print "\n\nOutput will be saved in the Output directory..."
	print "\n\nNote: These simulations take time, as each location requires 25050 ms of simulation data!"
	initsub()
	GenerateChirp()
	print "Chirp Generated....\n"
	count=0
	
	update_init()
	forsec Trunk {
		for(x) {
			if(x != 0) { // x=0 distance is equal to x=1 of prev section
				print count, " ", secname(), " ", x, raddist(x)
				st=new IClamp(x)
				st.dur=tstop
				st.del=0
				v3=new Vector()
				v3.record(&v(x))
			
				ChirpStim.play(&st.amp,dt)
			
				f1=new File()
			
				finitialize(v_init)
				fcurrent()
				while (t < tstop){
					fadvance()
				}
				
				sprint(str1,"Output/ChirpOutput_%d.txt",count)
				f1.wopen(str1)
				v3.printf(f1,"%f\n")
				f1.close()
				count=count+1
			}
		}
	}	
}	

/*****************************************************************/
// Calculate Input resistance of various trunk compartments.

objref st
objref f1

proc RN_Trunk () {
	print "\n\nComputing Input Resistances Along Trunk..."
	print "\n\nOutput will be saved in Trunk_RN.txt..."
	initsub()
	v3=new Vector()
	v3.record(&soma.v(0.5))
	f1=new File()
	tstop=250
	
	
	f1.wopen("Trunk_RN.txt")

	update_init()
	count=0
	forsec Trunk {
		for(x) {
			if(x != 0) { // x=0 distance is equal to x=1 of prev section

				st=new IClamp(0.5)
				st.dur=tstop
				st.del=0
				st.amp=-0.1
	
				v3=new Vector()
				v3.record(&v(x))
			
				finitialize(v_init)
				fcurrent()
	
				while (t < tstop){
					fadvance()
				}
				f1.printf("%f\t%f\n",raddist(x),-((v3.x[tstop/dt-1]-v_init)/100e-6)*1e-3)
				print secname(), " ", x," ", raddist(x)," ", -((v3.x[tstop/dt-1]-v_init)/100e-6)*1e-3
			}	
		}	
		count=count+1
	}	
	f1.close()
}

/*****************************************************************/


// Use the following code if you want to visualize the h channel gradient
// in a shape plot.

//objref sp
//sp=new PlotShape()
//sp.variable("ghdbar_hd")
//sp.show(0)
//sp.exec_menu("Shape Plot")
//sp.scale(1e-5,2e-4)

xpanel("Chirp Stimulus Responses in a 3D model with h Channels")
xlabel("Select one of the simulations below")
xbutton("Compute Input Resistances Along Trunk", "RN_Trunk()")
xbutton("Save Local Chirp Responses for Locations Along Trunk", "Chirp_Trunk()")
xpanel()


Loading data, please wait...