Active dendrites shape signaling microdomains in hippocampal neurons (Basak & Narayanan 2018)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:244848
The spatiotemporal spread of biochemical signals in neurons and other cells regulate signaling specificity, tuning of signal propagation, along with specificity and clustering of adaptive plasticity. Theoretical and experimental studies have demonstrated a critical role for cellular morphology and the topology of signaling networks in regulating this spread. In this study, we add a significantly complex dimension to this narrative by demonstrating that voltage-gated ion channels (A-type Potassium channels and T-type Calcium channels) on the plasma membrane could actively amplify or suppress the strength and spread of downstream signaling components. We employed a multiscale, multicompartmental, morphologically realistic, conductance-based model that accounted for the biophysics of electrical signaling and the biochemistry of calcium handling and downstream enzymatic signaling in a hippocampal pyramidal neuron. We chose the calcium – calmodulin – calcium/calmodulin-dependent protein kinase II (CaMKII) – protein phosphatase 1 (PP1) signaling pathway owing to its critical importance to several forms of neuronal plasticity, and employed physiologically relevant theta-burst stimulation (TBS) or theta-burst pairing (TBP) protocol to initiate a calcium microdomain through NMDAR activation at a synapse.
Reference:
1 . Basak R, Narayanan R (2018) Active dendrites regulate the spatiotemporal spread of signaling microdomains. PLoS Comput Biol 14:e1006485 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Dendrite; Synapse; Channel/Receptor; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampus CA1 pyramidal cell;
Channel(s): Ca pump; I A; I_SERCA; I Calcium; I_K,Na; I h; I Potassium;
Gap Junctions:
Receptor(s): AMPA; NMDA;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Active Dendrites; Detailed Neuronal Models; Calcium dynamics; Reaction-diffusion; Signaling pathways; Synaptic Plasticity;
Implementer(s): Basak, Reshma [reshmab at iisc.ac.in]; Narayanan, Rishikesh [rishi at iisc.ac.in];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal cell; AMPA; NMDA; I A; I h; I Calcium; I Potassium; I_SERCA; I_K,Na; Ca pump;
/
Basak_Narayayanan_2018
Spine1000_sample
readme.txt
apamp.mod
caminmax.mod
car.mod
cat.mod
ghkampa.mod
ghknmda.mod
h.mod
kadist.mod *
kaprox.mod
kdrca1.mod *
modcamechs.mod
na3.mod
nax.mod
vmax.mod
distance.hoc
distance_SD.hoc
Fig_13.hoc
mosinit.hoc
n123.hoc *
ObliquePath.hoc *
oblique-paths.hoc *
Sample_output_Calcium.txt
                            
//This HOC file was used to generate Figure 13 of the article. The parameters for generating other figures may be found in Table S1 of the article.
// The most important procedures in terms of adding calcium handling mechanisms, spines, recording species etc are:
// addspine()
//createspines()
//spineconnect()
//AddSynapse2()
//Record_species()
//defaultgrad_Dependent():run this to get the results in .txt format
//the calcium handling .mod file is called modcamechs, which is inserted in apical[59], as well as all the spines originating from it

//**************************************************************//
load_file("stdlib.hoc")
load_file("ObliquePath.hoc")
load_file("nrngui.hoc")
objref f2, f3, f4, f7
//**************************************************************//
// Current clamp Parameters
//**************************************************************//

objref cvode
cvode=new CVode()
cvode.active(0)

steps_per_ms = 40
dt = 0.025

stamp=2		// Stimulus Amplitude and Delay
stdur=1
stdel=2

//**************************************************************//
// Synaptic Parameters
//**************************************************************//
NAR=1.5
//P1 =20e-7//20e-7
//**************************************************************//
// Resistances and Capacitances
//**************************************************************//
//or just keep ra flat to 80 and gh just 20 fold increase
rasoma  = 120			
raend   = 70
raaxon = rasoma
ra=rasoma

rm     = 125000
rmsoma = rm
rmdend = rm
rmend=85000
rmaxon = rm

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

v_init    = -65
celsius   = 34
//**************************************************************//
// Active conductance densities and other related parameters
//**************************************************************//
//these values account for the bAP, Rin after getting ek=-90
Ek = -90
Ena = 55
Eh=-30
// for lovey16
gh=25e-6	// Base value of gh
caT=80e-6 
gna=0.016
gnais=gna	//0.008
gkdr=0.01
gka=0.0031
gkagrad=8
gka_pfactor=1
gka_dfactor=1
//**************************************************************//
// oblique conductances
//**************************************************************//
gka_oblique=0

//**************************************************************//
create axon[2]
objref sh, st, apc, vec, g[20], i_vec, v_vec, stdstt, fitvec
strdef str1, str2, str3
objref SD, AXON, SA, Basal, Trunk, AIS
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, st2
objref Trial
objref local_trunk
objref pl[150],opl[150]
objref v3, v4
strdef str1, rntfilename

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

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

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

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

proc set_Na_K_rad () {local xdist, ndist

	
	forsec SD {			
			insert na3 gbar_na3=gna ar2_na3=1
			insert kdr gkdrbar_kdr=gkdr
			insert kad gkabar_kad=gka
			insert kap gkabar_kap=gka
			for (x) {
				rdist=raddist(x)
				if(rdist<100){
						gkabar_kap(x)=gka*(1+gkagrad*rdist/100.0)
						gkabar_kad(x)=0
				} else {
					gkabar_kad(x)=gka*(1+gkagrad*rdist/100.0)
					gkabar_kap(x)=0
				}
			}
	}

	forsec "basal"{
			insert na3 gbar_na3=gna ar2_na3=1
			insert kdr gkdrbar_kdr=gkdr
			gkabar_kad=gka
			gkabar_kap=gka 
	}
  
	forsec AIS {			
			insert nax gbar_nax=gna*5
			insert kdr gkdrbar_kdr=gkdr
	}
    

	forsec "apical"{
		ar2_na3=0.8     // Add Na+ inactivation on apical side.
	}
	forall if (ismembrane("na3")||ismembrane("nax")) ena=Ena
	forall if (ismembrane("kap")||ismembrane("kad")||ismembrane("kdr")) ek=Ek
}

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

proc setactivess () {local xdist, ndist

	forsec SD {	// Trunk
		insert hd  

		for (x) {
			xdist=raddist(x)
			
			ghdbar_hd(x) = gh*(1+12/(1+exp(-(xdist-320)/50)))

			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
			}	

		}
	}
	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 set_cat() {local xdist, ndist

	forsec SD {	// Trunk
  			insert cat

		for (x) {
			xdist=raddist(x)
			gcatbar_cat(x)=caT*(1+30/(1+exp(-(xdist-350)/50)))


		}
	}
		forsec "basal"{
			insert cat gcatbar_cat=caT
	}
}

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

proc setactive() {
	setactivess ()
	set_Na_K_rad ()
	set_cat()
	add_vmax_mechanism()
	
	print "all active conductances set"
}	
/********************************************************************/

proc current_clamp() {
	access soma[0]
	st=new IClamp(0.5)
	st.amp=stamp
	st.del=stdel
	st.dur=stdur
}
/**********************************************************************/
objref ic[10]
proc current_clamp_train() {
    
	for (i=0; i<=9; i+=1){
		soma [0] ic[i] = new IClamp(0.5)
		ic[i].amp =stamp
		ic[i].del= 10 + i*25
		ic[i].dur= 2
	}
    
}
/**********************************************************************/

proc init_cell() {
	setpassive()
	addaxon()
    addspine()
    setactive()
	access soma[0]	// Reinitializing distance origin
	distance()
    
	finitialize(v_init)
	fcurrent()
        forall {
            for (x) {
			if (ismembrane("hd")||ismembrane("cat")||ismembrane("na3")||ismembrane("nax")||ismembrane("kdr")||ismembrane("kap")||ismembrane("kad")) {
		e_pas(x)=v(x)+(ica(x)+i_hd(x)+ina(x)+ik(x))/g_pas(x)


                } else {
        			e_pas(x)=v(x)
                }
            }
        }
}
/********************************************************************/

proc update_init(){
	finitialize(v_init)
	fcurrent()
        forall {
            for (x) {
				if (ismembrane("hd")||ismembrane("cat")||ismembrane("na3")||ismembrane("nax")||ismembrane("kdr")||ismembrane("kap")||ismembrane("kad")) {
				e_pas(x)=v(x)+(ica(x)+i_hd(x)+ina(x)+ik(x))/g_pas(x)


               } else {
        			e_pas(x)=v(x)
                }
            }
        }
}

/**********************************************************************/
create head1
create neck1
//creating the synapse-containing spine at the central location of apical[59] oblique and assign the channel types as well as the calcium handling mechanisms
proc addspine() {

    head1 {
        L = 0.5
        diam = 0.5
        nseg = 10
        insert pas
        insert na3
        insert kdr
        insert kad
        insert cat
        insert hd
        insert modcamechs
        
        cm = 1 // these values correspond to the oblique segment conductance value from which this spine originates
        g_pas = 8e-06
        Ra = 120
        e_pas=-65
        gbar_na3 = 0.16
        gkdrbar_kdr = 0.01
        gkabar_kad = 0.060547099
        gcatbar_cat = 0.00028570567
        ghdbar_hd = 6.8768407e-05
        
    }
    
    neck1 {
        L = 1
        diam = 0.1
        nseg = 20
        insert pas
        insert na3
        insert kdr
        insert kad
        insert cat
        insert hd
        insert modcamechs
        cm = 1
        g_pas = 8e-06
        Ra = 120
        e_pas=-65
        gbar_na3 = 0.16
        gkdrbar_kdr = 0.01
        gkabar_kad = 0.060547099
        gcatbar_cat = 0.00028570567
        ghdbar_hd = 6.8768407e-05
        
    }
    connect head1(1), neck1(0)
    connect neck1(1), apical[59](0.5) //connecting the synapse containing spine to the oblique of choice

}
/********************************************************************/
spineno=1000 //no. of spines on the oblique of choice (apical[59])
create spineheads[spineno]
create spinenecks[spineno]

//creating 1000 spines and connecting them to the oblique apical[59]
proc createspines(){
    for(a=0;a<spineno;a=a+1){
        spineheads[a] {
            print a
            L = 0.5
            diam = 0.5
            nseg = 1
            insert pas
            insert na3
            insert kdr
            insert kad
            insert cat
            insert hd
            insert modcamechs
        }
        spinenecks[a] {
            L = 1
            diam = 0.1
            nseg = 1
            insert pas
            insert na3
            insert kdr
            insert kad
            insert cat
            insert hd
            insert modcamechs
        }
        connect spineheads[a](1), spinenecks[a](0)
    }
    
}
/*****************************************************************/


proc load_3dcell() {

	forall delete_section()
	xopen($s1)

    create head1
    create neck1
    create spineheads[1000]
    create spinenecks[1000]
    
	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()

// here it is rad distacnes
//	soma[0] Trial.append()
//	apical[41] Trial.append() // 160 µm
	apical[71] Trial.append() // 320 µm
	apical[103] Trial.append() //417

	local_trunk = new SectionList()

	apical[60] local_trunk.append() 
	apical[62] local_trunk.append() 
	apical[64] local_trunk.append() 
	apical[65] local_trunk.append() 
	apical[69] local_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
	}
    
    apical[59]{
		nseg= 2000//2000
		update_init()
    }

	access soma[0]
	distance()

	init_cell()

}

/**********************************************************************/
// 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()
	AIS = new SectionList()

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

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

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


proc initsub() {
	load_3dcell("n123.hoc")
    handle_calmech ()
	finitialize(v_init)
	fcurrent()
}
/*****************************************************************/
 proc handle_calmech () {
 apical[59]{
  insert modcamechs
  insert caminmax
  }
 }
 
 /*****************************************************************/

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)
}

//**************************************************************//
// Alpha Function Parameters
//**************************************************************//

tau=10
nalpha=5
totalduration=500
alphaamp=0.1/exp(-1) // For normalization of the alpha amplitude to 1
alphafreq=20
alphadelay=50
/********************************************************************/
objref f1
proc AP_Amp () { local Comp
	f1=new File()
	//sprint(rntfilename,"Output/Linearity0_RN_%d.txt",gh*1e6)
	//f1.wopen(rntfilename)
	tstop=50
	count=0
	f1.wopen("1APAmp_Final.txt")
	current_clamp()

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

				finitialize(v_init)
				fcurrent()

				while (t < tstop){
				fadvance()
				}

				print count, " ", secname(), " ", x," ", raddist(x)," ", "\t", vm_vmax(x)-v_init
				f1.printf("%f\n", vm_vmax(x)-v_init)
				count=count+1
			}
		}
	}

	
	f1.close()
}
/**********************************************************************/
	
proc add_vmax_mechanism() {

	forall {
		insert vmax
	}
}
/**********************************************************************/

objref apamp[620]
objref apvec[620]
 
proc make_apamps() {local xval
	
	forall {
		insert vmax
	}

	ApAmpcnt=0
//	forsec Trunk{
	forsec Trial{
		for(x) {
//			if(x !=0){
			if(x ==0.5){
				apamp[ApAmpcnt]=new APAmp(x)
				//apvec[ApAmpcnt]=new Vector()
				ApAmpcnt=ApAmpcnt+1
			}
		}
	}
	
	//print ApAmpcnt

	for (i=0;i<ApAmpcnt;i+=1) {
		apamp[i].thresh=-64
		//apamp[i].record(apvec[i])
	}
}


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

objref f1
proc AP_Trunk () { local trunk
	f1=new File()
    tstop=50
	f1.wopen("Trunk_AP.txt")
	current_clamp()
	make_apamps()
	update_init()
			
	finitialize(v_init)
	fcurrent()
	for (i=0;i<ApAmpcnt;i+=1) {
		apamp[i].max=v_init
	}
	
	while (t < tstop){
		fadvance()
	}

	for (i=0;i<ApAmpcnt;i+=1) {
		if (apamp[i].max>=-64){
                	f1.printf("%f\n", apamp[i].max-v_init)
			print i, " ", apamp[i].max-v_init

        	} else {
			print "Propagation failure!!!!  ", i
		}
	}
	
	f1.close()
}

/*****************************************************************/
objref st
objref f1

proc RN_Trunk_Trial () {
	f1=new File()
	tstop=500
	cvode.active(1)
	f1.wopen("Trunk_RN_basemod_trial.txt")
	update_init()
	count=0
	forsec Trial {
		for(x) {
			if(x == 0.5) { // x=0 distance is equal to x=1 of prev section

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

/*****************************************************************/
objref st
objref f1

proc RN_Trunk () {
	f1=new File()
	tstop=300
	cvode.active(1)
	
	f1.wopen("Trunk_RN_basemod.txt")
	

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

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



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

proc RN_Fit () {local trunk
	f1=new File()
	i_vec=new Vector(11) //current vector
	stdstt=new Vector(11)	//stdstt=steady state
	fitvec=new Vector(11)	//for fitting
	cvode.active(1)

	tstop=500
			
	f1.wopen("2Rin_CaT_Final.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(x)
				for(i=0;i<=10; i+=1){
					st.dur=tstop
					st.del=50
					st.amp=((i-5)*10)/1000
					i_vec.x[i]=(i-5)*10/1000
					v_vec=new Vector(11) //voltage
					v_vec.record(&v(x))

			
					finitialize(v_init)
					fcurrent()
	
					while (t < tstop){
						fadvance()
					}
					stdstt.x[i]=(v_vec.x[v_vec.size()-1]-v_init)
				}	
				a=stdstt.x[0]/i_vec.x[0]
				b=0

				error=stdstt.fit(fitvec, "line", i_vec, &a, &b)
				print secname(), " ", x," ", raddist(x)," ", a
				f1.printf("%f\n",a)
				count=count+1
				
			}
		}
	}
f1.close()

}
/********************************************************************/

proc Recorcd_bAP_Vm () {local trunk
	tstop=50
	update_init()
	count=0
	current_clamp()
	apical[71] {

				v3=new Vector()
				v3.record(&v(0.5))
					
				f1=new File()
			
				finitialize(v_init)
				fcurrent()
				while (t < tstop){
					fadvance()
				}
				
				sprint(rntfilename,"Output/Optimized_Cell2/bAP_Trunk/TrunkbAP_%d.txt",300)

				f1.wopen(rntfilename)
				v3.printf(f1,"%f\n")
				f1.close()
				count=count+1
		}
}	
/********************************************************************/

proc Recorcd_Vm () {local trunk
	tstop=550
	update_init()
	count=0
	apical[71] {
		
				st=new IClamp(0.5)
				st.dur=300
				st.del=50
				st.amp=-0.05

				v3=new Vector()
				v3.record(&v(x))
					
				f1=new File()
			
				finitialize(v_init)
				fcurrent()
				while (t < tstop){
					fadvance()
				}
				
				sprint(rntfilename,"Output/Optimized_Cell2/RN_Trunk/TrunkVmH_%d.txt",300)

				f1.wopen(rntfilename)
				v3.printf(f1,"%f\n")
				f1.close()
				count=count+1
		}
}	
/********************************************************************/

proc Recorcd_Vm_RN () {local trunk
	tstop=550
	update_init()
	count=0

	forsec Trial {
		for(x) {
			if(x == 0.5) { // x=0 distance is equal to x=1 of prev section
				st=new IClamp(x)
				for(i=-50;i<=50; i+=10){
					st.dur=300
					st.del=50
					st.amp=i/1000

					v3=new Vector()
					v3.record(&v(x))
						
					f1=new File()
				
					finitialize(v_init)
					fcurrent()
					while (t < tstop){
						fadvance()
					}
				
				sprint(rntfilename,"Output/Vm_Trace_RN/Loc_Amp_%d_%d.txt",count,i)

				f1.wopen(rntfilename)
				v3.printf(f1,"%f\n")
				f1.close()
				}
				count=count+1
			}	
		}
	}	
}	

/********************************************************************/
objectvar Ecellvec, Ecellvec1, f2, FRVec
strdef rntfilename1, rntfilename2
proc Firing_Rate () {local trunk
	tstop = 1100
	soma[0]{
		apc=new APCount(0.5)
		apc.thresh=15
		Ecellvec = new Vector()
		v3 = new Vector()
		FRVec = new Vector(6)
		count=0
		
	for(i=0;i<=250; i+=50){
				f1=new File()
				f2=new File()
				st=new IClamp(0.5)
				st.dur=1000	//tstop
				st.del=100
				st.amp=i/1000
							
				finitialize(v_init)
				fcurrent()
		
				while (t < tstop){
					fadvance()
				}
				FRVec.x[count]=apc.n
				print apc.n
			
				
			count+=1
		}
				sprint(rntfilename1,"Output/Firing_Rate.txt")
				f1.wopen(rntfilename1)
				FRVec.printf(f1,"%f\n")
				
				f1.close()
		
		}		
}


/*****************************************************************/
proc get_gh(){
	cc=0
	for(yy=10; yy<=60; yy+=10){
		print yy
		set_h()
		cc+=1
		update_init()
		ZAP_Single()
	}

}
/*****************************************************************/
objref randspine
spinenumber = 1000 // no. of spines
segnumber = 2002
double spineloc [segnumber]
double xlist [segnumber]
double xindex[spinenumber]
double spine_cap[spinenumber]
double spine_memres[spinenumber]
double spine_axres[spinenumber]
double spine_restmem[spinenumber]
double spine_na3[spinenumber]
double spine_kdr[spinenumber]
double spine_ka[spinenumber]
double spine_cat[spinenumber]
double spine_h[spinenumber]

// This procedure connects 1000 spines in a randomly distributed fashion in the chosen oblique (apical[59]) and assigns all membrane properties of the oblique segment to each each spine originating from that segment

proc spineconnect(){
    createspines()
    for(i=0;i<segnumber;i+=1){
        spineloc[i]=0
    }
    randspine = new Random(11)
    randspine.uniform(0,segnumber)
    
    count = 0
    
    while(count<spinenumber){
        x1 = int(randspine.repick())
        while(spineloc[x1]==1){
            x1 = int(randspine.repick())
        }
        spineloc[x1]=1
        count+=1
    }
    
    count=0
    apical[59]{
        for(x){
            xlist[count]=x
            count+=1
        }
    }
    count = 0
    for(i=0;i<segnumber;i+=1){
        if (spineloc[i]==1){
            xindex[count]=xlist[i]
            count+=1
        }
    }
    access apical[59]
    count1=0
    apical[59]{
        for(x){
            if(x==(xindex[count1])){
                connect spinenecks[count1](1), apical[59](x)
                spine_cap[count1]= 1
                spine_memres[count1]= g_pas(x)
                spine_axres[count1]= Ra
                spine_restmem[count1]= e_pas(x)
                spine_na3[count1]= gbar_na3(x)
                spine_kdr[count1]= gkdrbar_kdr(x)
                spine_ka[count1]= gkabar_kad(x)
                spine_cat[count1]= gcatbar_cat(x)
                spine_h[count1]= ghdbar_hd(x)
                if (count1<spinenumber-1){
                    count1+=1
                }
            }
        }
    }
    
    for(i=0;i<spinenumber;i+=1){
        spineheads[i]{
            for(x){
                cm = spine_cap[i]
                g_pas = spine_memres[i]
                Ra = spine_axres[i]
                e_pas = spine_restmem[i]
                gbar_na3 = spine_na3[i]
                gkdrbar_kdr = spine_kdr[i]
                gkabar_kad = spine_ka[i]
                gcatbar_cat = spine_cat[i]
                ghdbar_hd = spine_h[i]
            }
        }
        spinenecks[i]{
            for(x){
                cm = spine_cap[i]
                g_pas = spine_memres[i]
                Ra = spine_axres[i]
                e_pas = spine_restmem[i]
                gbar_na3 = spine_na3[i]
                gkdrbar_kdr = spine_kdr[i]
                gkabar_kad = spine_ka[i]
                gcatbar_cat = spine_cat[i]
                ghdbar_hd = spine_h[i]
            }
        }
    }
}


/*****************************************************************/
objref f1, f2, netlist, s[1], f1, DEND, sapamp, somavec, sampvec, f2
strdef  str2
objref s, ampasyn, nmdasyn
strdef rntfilename2
//this procedure records the voltage generated at the soma with input to each synapse. Useful for determining what permeability value to use for each synapse
proc find_epsp_amplitudes() { local x
	netlist=new List()
    
	tstop=40
    Gstart=3e-8
    Gend=1
    f1=new File()
    f2=new File()
	sprint(rntfilename,"EPSPAmps_synapse_spine.txt")
    sprint(rntfilename2,"EPSP_synapse_spine.txt")
	f1.wopen(rntfilename)
    f2.wopen(rntfilename2)
	soma[2] s = new NetStim(0.5)
	s.interval=1   // ms (mean) time between spikes
	s.number=1     // (average) number of spikes
	s.start=2   // ms (mean) start time of first spike
	s.noise=0      // range 0 to 1. Fractional randomness.
    
    countsegment =0
    head1 {
            for (x) {
                if(x!=0 && x<1){
                    
                    G=Gstart
                    ampasyn=new ghkampa(x)
                    ampasyn.taur=2
                    ampasyn.taud=10
                    ampasyn.Pmax=Gstart
                    netlist.append(new NetCon(s,ampasyn,0,0,0))

                    while(G<=Gend){
                        ampasyn.Pmax=G
                        update_init()
                        
                        MAX=-65
                        while (t <tstop){
                            fadvance()
                            //print soma[2].v(0.5)
                            if(soma[2].v(0.5)>MAX) {
                                MAX=soma[2].v(0.5)
                            }
                        }
                        print secname(), "\t ", x, "\t", " ", G, " ", MAX, "\t" , "\t", 65+MAX
                        countsegment=countsegment+1
                        if(MAX>=-64.8) {
                            print secname(),"reached", " ", x, " ", " ", G, " ", MAX, "\t" , "\t", 65+MAX
                            f1.printf("%g\t%g\t%g\n",x,G, MAX)
                            f2.printf("%g\n", G)
                            update_init()
                            break
                        }
                        G=G+5e-7
                    }
                    if(G>=Gend){
                        print secname(), " ", x," ", -1,  " ", G, " ", MAX, "\t" , "\t", 65+MAX
                    }
                    
                }
    
        }
    }
    print countsegment
    f1.close()
    f2.close()
}
/*****************************************************************/
//added_reshma
objref ampa, nmda, spikegen, ncl
objref calcvec[1005],camvec[1005],camca4vec[1005],camkiivec[1005],camkiicamca4vec[1005],pcamkiivec[1005],PP1vec[1005],dPP1vec[1005],voltagevec[1005]

objref f, f1, f5, f6, f8, f9, f10, f13, f14
strdef filename, filename1, filename3, filename4, filename5, filename6, filename7, filename10, filename11
create dend, oblique1, apical[59], oblique3
objref ampa2[1],nmda2[1], ncl3[1], ns2[150]
/*************************************************************************/

//this procedure adds the NMDA and AMPA synapse at the spine head located at the middle of the chosen oblique:apical[59]
proc AddSynapse2() {
    
    tstop=10000
    
    for (i=0;i<30;i+=1){
        //time=(20+i*200)
        //print "time=", time
        head1 ns2[i]= new NetStim(0.3)
        
        ns2[i].interval=10 // ms (mean) time between spikes
        ns2[i].number=5    // (average) number of spikes
        ns2[i].start= (20+i*200)
        ns2[i].noise=0      // range 0 to 1. Fractional randomness.
        
    }
    
    for (i=0;i<=0;i+=1){
        ampa2[i]= new ghkampa()         // create synapse
        head1 ampa2[i].loc(0.5)
        ampa2[i].taur = 2
        ampa2[i].taud = 10
        ampa2[i].Pmax = 1.53e-06
        nmda2[i] = new ghknmda ()
        head1 nmda2[i].loc(0.5)
        nmda2[i].taur = 5
        nmda2[i].taud = 50
        nmda2[i].Pmax = 1.53e-06 * NAR//NAR stands for NMDA:AMPA ratio
        
    }
    
    
    head1 ncl3[0]=new List()
    for (i=0;i<30;i+=1){
		head1 ncl3[0].append(new NetCon(ns2[i], ampa2[0], 0, 0, 0.0001))
        head1 ncl3[0].append(new NetCon(ns2[i], nmda2[0], 0, 0, 0.0001))
    }
}
    /********************************************************************************/
    
//This procedure creates the vectors necessary for storing the voltage, calcium, as well as the downstream molecules

objref temp_ca[1005], temp_cam[1005],temp_camca4[1005],temp_camkii[1005],temp_camkiicamca4[1005],temp_pcamkii[1005],temp_dPP1[1005],temp_PP1[1005]

proc Record_species(){ local i
    
    segcount=2000
    print "segcount=", segcount
    
    veccount=0
    i=0
    tcount=1
    FIRST=0
    
    access apical[59]
    for(i=0; i<1005; i+=1){         // creating temporary vectors to store all 40 pnts per millisecond for all the species except voltage
        temp_ca[i]=new Vector(40)
        temp_cam[i]=new Vector(40)
        temp_camca4[i]=new Vector(40)
        temp_camkii[i]=new Vector(40)
        temp_camkiicamca4[i]=new Vector(40)
        temp_pcamkii[i]=new Vector(40)
        temp_dPP1[i]=new Vector(40)
        temp_PP1[i]=new Vector(40)
    }
    
    i=0
    
    
	//for (x=1000/segcount; x<=1500/segcount; x+=1/segcount) { //while recording voltage, use this, segment 500 to 1000, then 1000 to 1500, otherwise memory allocation error my occur
    for (x=500/segcount; x<=1500/segcount; x+=1/segcount) { //while recording species other than voltage, use this, for all the 1000 segments around the synptic location

        print x
        calcvec[i]= new Vector() // creating vectors for saving all species
        camvec[i]= new Vector()
        camca4vec[i]= new Vector()
        camkiivec[i]= new Vector()
        camkiicamca4vec[i]= new Vector()
        pcamkiivec[i]= new Vector()
        PP1vec[i]= new Vector()
        dPP1vec[i]= new Vector()
//        voltagevec[i]= new Vector()
//        voltagevec[i].record (&apical[59].v(x)) // 40 voltage values are recorded every millisecond, to account for faster voltage fluctuations compared to calcium and other downstream species

        i=i+1
        veccount+=1
    }
    
    print "Veccount=", veccount
}
/********************************************************************/

// This function/procedure calls all the previous procedures relevant for handling the calcium mechanisms, adding the spine containing synapse,adding spines, recording all species and finally stores them in .txt files

proc defaultgrad_Dependent(){
    //calling all necessary functions
    AddSynapse2()
    spineconnect()
    Record_species()
    tstop= 100//10000 // total run time of simulation
    
    //Creating file handles for storing species values
    f=new File()
    f1=new File()
    f5=new File()
    f6=new File()
    f8=new File()
    f9=new File()
    f10=new File()
    f13=new File()
    f14=new File()
    
    
    reccount=0
    //mscount=0
    update_init()
    //print "Starting ka, Cat,h: ", gkabar_kad,gcatbar_cat,ghdbar_hd
    
        
  while (t < tstop) {
    if (reccount==39) {
        for (i=0;i<1000;i+=1){                          // storing mean of all 40 values of each species generated every millisecond
            calcvec[i].append(temp_ca[i].mean())
            camvec[i].append(temp_cam[i].mean())
            camca4vec[i].append(temp_camca4[i].mean())
            camkiivec[i].append(temp_camkii[i].mean())
            camkiicamca4vec[i].append(temp_camkiicamca4[i].mean())
            pcamkiivec[i].append(temp_pcamkii[i].mean())
            PP1vec[i].append(temp_PP1[i].mean())
            dPP1vec[i].append(temp_dPP1[i].mean())
        }
        reccount=0
        //mscount+=1
        //print mscount, calcvec[500].x[mscount-1]
    }
   
    i=0
    apical [59] {
        for (x) {
            if (i>=500 && i<1500){
                temp_ca[i-500].x[reccount]=apical[59].cai(x)        //storing 40 pnts per millisecond for all the species except voltage in the temporary vectors created above
                temp_cam[i-500].x[reccount]=apical[59].cam_modcamechs[0](x)
                temp_camca4[i-500].x[reccount]=apical[59].camca4_modcamechs[0](x)
                temp_camkii[i-500].x[reccount]=apical[59].camkii_modcamechs[0](x)
                temp_camkiicamca4[i-500].x[reccount]=apical[59].camkii_camca4_modcamechs[0](x)
                temp_pcamkii[i-500].x[reccount]=apical[59].pcamkii_camca4_modcamechs(x)
                temp_PP1[i-500].x[reccount]=apical[59].PP1_modcamechs(x)
                temp_dPP1[i-500].x[reccount]=apical[59].dPP1_modcamechs(x)
                
            }
            i=i+1
        }
    }
    reccount+=1
    //print reccount
    fadvance()
}

print "Done Running "
    
         for (i=0; i<veccount; i+=1) { // saving all the vectors to files
            
            sprint (filename,"Spine1000_sample/Calcium_%d.txt",i)
            f.wopen(filename)
            calcvec[i].printf(f,"%g\n")
            f.close()
        
            sprint (filename4,"Spine1000_sample/Cam_%d.txt",i)
            f6.wopen(filename4)
            camvec[i].printf(f6,"%g\n")
            f6.close()
        
            sprint (filename1,"Spine1000_sample/Camca4_%d.txt",i)
            f1.wopen(filename1)
            camca4vec[i].printf(f1,"%g\n")
            f1.close()
        
            sprint (filename5,"Spine1000_sample/Camkii_%d.txt",i)
            f8.wopen(filename5)
            camkiivec[i].printf(f8,"%g\n")
            f8.close()
        
            sprint (filename3,"Spine1000_sample/CamCamkii_%d.txt",i)
            f5.wopen(filename3)
            camkiicamca4vec[i].printf(f5,"%g\n")
            f5.close()
        
            sprint (filename10,"Spine1000_sample/pCamkii_%d.txt",i)
            f13.wopen(filename10)
            pcamkiivec[i].printf(f13,"%g\n")
            f13.close()
        
            sprint (filename6,"Spine1000_sample/PP1_%d.txt",i)
            f9.wopen(filename6)
            PP1vec[i].printf(f9,"%g\n")
            f9.close()
             
            sprint (filename7,"Spine1000_sample/dPP1_%d.txt",i)
            f10.wopen(filename7)
            dPP1vec[i].printf(f10,"%g\n")
            f10.close()
        
//            sprint (filename11,"Spine1000_sample/Voltage_%d.txt",i+500)//change this when u change voltage rec location0->500
//            f14.wopen(filename11)
//            voltagevec[i].printf(f14,"%g\n")
//            f14.close()
        }
    print "Done Saving"
	
}
/*****************************************************************/

//calling the necessary functions
initsub()
defaultgrad_Dependent()












Loading data, please wait...