High frequency oscillations in a hippocampal computational model (Stacey et al. 2009)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:135902
"... Using a physiological computer model of hippocampus, we investigate random synaptic activity (noise) as a potential initiator of HFOs (high-frequency oscillations). We explore parameters necessary to produce these oscillations and quantify the response using the tools of stochastic resonance (SR) and coherence resonance (CR). ... Our results show that, under normal coupling conditions, synaptic noise was able to produce gamma (30–100 Hz) frequency oscillations. Synaptic noise generated HFOs in the ripple range (100–200 Hz) when the network had parameters similar to pathological findings in epilepsy: increased gap junctions or recurrent synaptic connections, loss of inhibitory interneurons such as basket cells, and increased synaptic noise. ... We propose that increased synaptic noise and physiological coupling mechanisms are sufficient to generate gamma oscillations and that pathologic changes in noise and coupling similar to those in epilepsy can produce abnormal ripples."
Reference:
1 . Stacey WC, Lazarewicz MT, Litt B (2009) Synaptic noise and physiological coupling generate high-frequency oscillations in a hippocampal computational model. J Neurophysiol 102:2342-57 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampus CA1 pyramidal GLU cell; Hippocampus CA3 pyramidal GLU cell; Hippocampus CA1 interneuron oriens alveus GABA cell; Hippocampus CA1 basket cell;
Channel(s): I Na,t; I A; I K; I h;
Gap Junctions: Gap junctions;
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Oscillations;
Implementer(s): Lazarewicz, Maciej [mlazarew at gmu.edu]; Stacey, William [wstacey at med.umich.edu];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal GLU cell; Hippocampus CA3 pyramidal GLU cell; Hippocampus CA1 interneuron oriens alveus GABA cell; GabaA; AMPA; NMDA; I Na,t; I A; I K; I h;
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE
//
// Copyright 2010, The University Of Michigan
// 	
//   All rights reserved.
//   For research use only; commercial use prohibited.
//   No Distribution without permission of William Stacey
//   wstacey@umich.edu
//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

begintemplate TGnet

external cvode

// PUBLIC VARIABLES

// PUBLIC REFERENCES
public pnm   

// PUBLIC METHODS
public setScatteredVoltages, activeSynapses, pr, signalSynapse, activeSynapsesRandom, writeVoltages, recordVoltages, noiseStim, nc, nil, noiselist, signallist, par_gaps //, ncstim
public randomPyrThenStim, writesingleVoltage, writenoise, recordnoise, manualParams

objref pnm  
objref synParamSet, synCellSet, noisepyrlist, noisebasklist, nc, nil, noiselist, signallist, par_gaps

	manParamFlag=0
	gmax_basktopyr=0
	gmax_pyrgap=0   //For manual parameter entry


{load_file("PYRkop.tem")}
{load_file("Bwb.tem")}
{load_file("Ok.tem")}
//{load_file("NoiseStim.tem")}
{load_file("../parameters/synapses.tem")}
{load_file("../parameters/cells.tem")}



// =================================================================================================
//
// init(Npyr, Nbask, Nolm)
//
// =================================================================================================
proc init() {
	
	print manParamFlag, "entered init"

	oneCell = 0   //I don't know 
	nModules  = 1    // set up for # of cell columns, was 4 in paper  ###changed to 1 on 3/11/08
	debugFlag = 0
	
	synCellSet = new CellParamSet()  //reads in the cells.par into synCellSet
	
    Npyr  = synCellSet.cellSet.object(0).N + 1
    Nbask = synCellSet.cellSet.object(1).N
    Nolm  = synCellSet.cellSet.object(2).N
    Nnoise  = Npyr + Nbask //one noise synapse on each pyr and bask
    
    Npyr  = Npyr      
    Nbask = Nbask 
    Nolm  = Nolm  
    Nnoise = Nnoise  //added
    
    N     = Npyr+Nbask+Nolm//+Nnoise   //added the Nnoise
   
	pnm = new ParallelNetManager(nModules*N)
   
    if (debugFlag) print "Assigning GIDs"
	//assignGIDs()
	pnm.round_robin()
	
	if (debugFlag) print "Creating cells"
    createCells()
      if (debugFlag) print "Connecting network"

    connectNetwork()
    if(nModules>1) connectNetworkInterModules()
    //activeSynapses()
	activeSynapsesZero()
	signallist=new List()
	
}




// =================================================================================================
//
// iterators
//		pyrIter
//		baskIter
//		olmIter
//
// =================================================================================================
iterator pyrIter() { local i, ii

	for ii = 0, nModules-1 {
		 
		for i=0, Npyr-1 {
		
			$&1 = i
			$&2 = i+ii*N
			iterator_statement
		}
	}
//if (debugFlag) print "pyrIter"
}

iterator baskIter() { local i, ii

	for ii = 0, nModules-1{
	for i=0,Nbask-1 {
		
		$&1 = i
		$&2 = i+Npyr+ii*N
		iterator_statement
	}
	}
//if (debugFlag) print "baskIter"
}

iterator olmIter() { local i, ii

	for ii = 0, nModules-1{
	for i=0,Nolm-1 {
		
		$&1 = i
		$&2 = i+Npyr+Nbask+ii*N
		iterator_statement
	}
	}
//if (debugFlag) print "OLMiter"
}

iterator noiseIter() { local i, ii

	for ii = 0, nModules-1{
	for i=0,Nnoise-1 {
		
		$&1 = i
		$&2 = i+Npyr+Nbask+Nolm+ii*N
		iterator_statement
	}
	}
//if (debugFlag) print "Noiseiter"
}



// =================================================================================================
//
// assignGIDs()
// not called, rather the round_robin is called
// =================================================================================================
proc assignGIDs() { local i, gid
	
	i   = 0
	gid = 0
	
	if (debugFlag) print "ID#", pnm.pc.id, " NHOST:", pnm.pc.nhost
	
	if (pnm.pc.nhost>1) {
	
		for pyrIter(&i, &gid)   pnm.set_gid2node(gid, i%(pnm.pc.nhost-1))
		for baskIter(&i,&gid)   pnm.set_gid2node(gid, pnm.pc.nhost-1)
		for olmIter(&i, &gid)   pnm.set_gid2node(gid, pnm.pc.nhost-1)
		for noiseIter(&i, &gid)   pnm.set_gid2node(gid, pnm.pc.nhost-1)

		
	}else{
		
		for pyrIter(&i, &gid)   pnm.set_gid2node(gid, 0)
		for baskIter(&i,&gid)   pnm.set_gid2node(gid, 0)
		for olmIter(&i, &gid)   pnm.set_gid2node(gid, 0)
		for noiseIter(&i, &gid) pnm.set_gid2node(gid, 0)
	}
}



// =================================================================================================
//
// createCells()
//
// =================================================================================================
proc createCells() {local gc, i, pNN localobj parRef

	synParamSet = new SynParamSet(pnm.pc.id)	
		pNN = 0			
		for i=0,  synCellSet.cellSet.count()-1 {				
			if(strcmp(synCellSet.cellSet.object(i).name, "Pyr")==0) {							
				parRef = synCellSet.cellSet.object(i)
				pNN = pNN + 1
			}			
		}				
		if(pNN!=1) print "ERROR pNN PYR"	
	createPyrCells(parRef)
		pNN = 0		
		for i=0,  synCellSet.cellSet.count()-1 {				
			if(strcmp(synCellSet.cellSet.object(i).name, "Bask")==0) {							
				parRef = synCellSet.cellSet.object(i)
				pNN = pNN + 1
			}			
		}				
		if(pNN!=1) print "ERROR pNN BASK"	
	createBaskCells(parRef)	
			pNN = 0			
		for i=0,  synCellSet.cellSet.count()-1 {				
			if(strcmp(synCellSet.cellSet.object(i).name, "OLM")==0) {							
				parRef = synCellSet.cellSet.object(i)
				pNN = pNN + 1
			}			
		}				
		if(pNN!=1) print "ERROR pNN OLM"	
	createOLMCells(parRef)
}


// =================================================================================================
//
// createPyrCells( CellParam (object))
//
// =================================================================================================
proc createPyrCells() {local i, j, gid localobj c, r
	r = new Random(pnm.pc.id+startsw())	
	i   = 0
	gid = 0
    for pyrIter(&i, &gid) {    	
    	if (pnm.gid_exists(gid)) {		
		if(gid%N==0) { // Passive cell			
				//c = new PYRkop(gid, 0, $o1.iappSsd, 0, $o1.iappAsd, $o1.iappUnits)
				c = new PYRkop(gid, 0, 0, 0, 0, $o1.iappUnits)
				
//### added stim to first 20 cells 7/8/8
		//}else if(gid<20){
				//c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform($o1.iappAl, $o1.iappAh), $o1.iappAsd, $o1.iappUnits)
				//c = new PYRkop(gid, 0, 94.2, 180, 94.2, 0)				
		}else{			
			c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform($o1.iappAl, $o1.iappAh), $o1.iappAsd, $o1.iappUnits)
			if (0) {
				if (gid%N<12/40*N && int(gid/N)<2 ) {					
					c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform(330, 349), $o1.iappAsd, $o1.iappUnits)
					//c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform(300, 310), $o1.iappAsd, $o1.iappUnits)
				
				} else {				
					if (gid%N<12/40*N && int(gid/N)>=2) {					
						//c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform(820, 849), $o1.iappAsd, $o1.iappUnits)			
						c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform(350, 400), $o1.iappAsd, $o1.iappUnits)			
					
					}else{					
						c = new PYRkop(gid, r.uniform($o1.iappSl, $o1.iappSh)+0*(gid*80/400+250)*2.5e-1, $o1.iappSsd, r.uniform($o1.iappAl, $o1.iappAh), $o1.iappAsd, $o1.iappUnits)
					}
				}
			}
		}
		
		pnm.register_cell(gid, c )
		addSyn(synParamSet, gid, "Pyr")		
    }
}
}


// =================================================================================================
//
// createBaskCells(cell param object)
//
// =================================================================================================
proc createBaskCells() {local i, j, gid localobj c, r
	
	r = new Random(pnm.pc.id+startsw())
	
	i      = 0
	gid    = 0

    for baskIter(&i, &gid) if (pnm.gid_exists(gid)) {
		
		c = new Bwb(gid, r.uniform($o1.iappSl, $o1.iappSh), $o1.iappSsd, $o1.iappUnits)
	
		pnm.register_cell(gid, c )

	 	//if (debugFlag) print "Adding Basket synapses"
		addSyn(synParamSet, gid, "Bask")
		
		c.recordVoltage()
    }
}




// =================================================================================================
//
// createOLMCells(men, sd)
//
// =================================================================================================
proc createOLMCells() {local i, j, gid localobj c, r
	r = new Random(pnm.pc.id+startsw())
	i   = 0
	gid = 0

    for olmIter(&i, &gid) if (pnm.gid_exists(gid)) {
		
		//c = new Ok(gid,  r.uniform($o1.iappSl, $o1.iappSh), $o1.iappSsd, $o1.iappUnits)
		c= new Ok(gid, 0,0,0)
		pnm.register_cell(gid, c )

		addSyn(synParamSet, gid, "OLM")
		if (debugFlag) "creating OLM cells"
    }
}

// =================================================================================================
//
// createNoiseStim()
//
// =================================================================================================
proc createNoiseStim() {local i, j, gid localobj c 
	i   = 0
	gid = 0

    for noiseIter(&i, &gid) if (pnm.gid_exists(gid)) {
		
		c = new NetStim()

		pnm.register_cell(gid, c )
// as of 5/26/08, the noise cells are 121-201 (to pyr) and 202-221 (to bask)
//as of 6/08 this is no longer used, and the netstims are just NULLOBJECTs	
    }
}



iterator ltr() {local i, cnt
	for i = 0, $o2.count - 1 {
		$o1 = $o2.object(i)
		iterator_statement
	}
}



strdef cmd
// =================================================================================================
//
// activeSynapsesZero()
//
// =================================================================================================
proc activeSynapsesZero() {localobj xo, co
	
	co = pnm.nclist
	
	//if ($1) print "synapse ACTIVATION at ", pnm.pc.id, co.count()
	for ltr(xo, co) {
		
		//if ($1) print pnm.pc.id, xo.precell, xo.postcell, xo.syn, xo.weight
		sprint(cmd, "%s", xo.postcell)
		
		if(strcmp(cmd,"NULLobject")!=0) xo.active(0)
	}
}
objref synR
// =================================================================================================
//
// randomPyrThenStim()   to have the pyr cells firing randomly and lightly up to certain time, 
//      then start some input at that time
// =================================================================================================


proc randomPyrThenStim() {localobj xo, co, noiseR, noiseA, firstR
// $1 start stim, $2 end stim, $3 min pyrthresh (uniform) prior to stim
// $4 '0' DC current, '1' random inputs,  $5 current value/pyrthresh value

	noiseR = new Random( startsw() )
	synR = new Random( startsw() )
	firstR = new Random( startsw() )
	noiseA = new Random( startsw() )
	synR.uniform(5,15)
	co   = pnm.nclist
	for ltr(xo, co) {
		//print pnm.pc.id, xo.precell, xo.postcell, xo.syn, xo.weight, xo
		sprint(cmd, "%s", xo.postcell)
		if(strcmp(cmd, "NULLobject")!=0) {
			sprint( cmd, "%s.active(1)", xo )
			cvode.event(synR.repick(), cmd)    //makes NetCon[#] active at random times within 250, 750
			//###but does that make an event?? I don't think so!  
			// but if you remove this, there are no conns to the baskets...
		}
	}
//add a noise event to every pyr cell  and every bask cell
noiseR.uniform(0,100)
noiseA.poisson(0.8)
minthreshold=$3  //maximum pyrthresh on a uniform scale during first segment
// no basketthreshold is implemented now
basketthreshold=100
noisestepper=.1  //how long the time intervals are in ms between each poisson noise iteration
rng=0
tStop=$2  //this is where the entire run ends
tStart=$1  //this is where the stim starts, and where the noise ends
randomflag=$4 //1 will go to random inputs, 0 will be DC input

noiselist=new List()
	noisetime=0

// set up the Netcons for noise synapses on pyramidal cells
for ii=0, nModules-1 {
   for i=0,80 {
	noisetime=0
	gidnoise= i+ ii*N
    if(pnm.gid_exists(gidnoise)) {
		nc=new NetCon(nil, pnm.pc.gid2obj(gidnoise).synlist.o(2))
	//print gidnoise, pnm.pc.gid2obj(gidnoise).synlist.o(2), nc.postcell, nc
	nc.weight=.00001 * 53.407075   //cell area of pyr vs. the bask
	noiselist.append(nc)
	coweight=nc.weight
	firstR.uniform(minthreshold,100)
	threshold=firstR.repick()
	//print threshold, gidnoise
	while(noisetime<tStop){
		if (noisetime>0 && noisetime<tStart) {   //first segment, just noise
		  rng=noiseR.repick()
		  if (rng>threshold) {
			amp=noiseA.repick()
			noiselist.o[i].weight = nc.weight * amp
			//print noiselist.o[i], noisetime, noiselist.o[i].weight 
			sprint(cmd, "%s.active(1)",noiselist.o[i])
			cvode.event(noisetime,cmd)
			//print i, noisetime
			noiselist.o[i].event(noisetime)
		    } //rng	
		    nc.weight=coweight  //reset weight to original value, else it retains it from this event
		  } // if noisetime
		if (noisetime>tStart && noisetime<tStop) {  //second segment
			if (!randomflag) {
				//pnm.pc.gid2obj(gidnoise).setIamp(0,0,$5,0,0)  //( iB, iS, iA1, iA2, iA3 [uA/cm2])
				//that one goes constantly with no delay
				//pnm.pc.gid2obj(gidnoise).iappA1.set($5,tStart)  //could work but doesn't have area info
				pnm.pc.gid2obj(gidnoise).setAdend1Iamp($5,tStart) //40 saturates it, 20 is poorly coordinated
					//30 has a couple of saturated ones, 35 a few more
			} else {
				rng=noiseR.repick()
		  		if (rng>$5) {  //now all receive same input, as def by $5=pyrthresh
					amp=noiseA.repick()
					noiselist.o[i].weight = nc.weight * amp
					sprint(cmd, "%s.active(1)",noiselist.o[i])
					cvode.event(noisetime,cmd)
					noiselist.o[i].event(noisetime)
		    		} //rng	
		    		nc.weight=coweight  
			} //randomflag
		} //if noisetime
		

		noisetime=noisetime+noisestepper
		} //while noisetime
    } //gidexists
	} // for i
//print noiselist.count
//print noiselist.o[1].postcell, noiselist.o[1].weight, noiselist.o[1]
   for i=81,100 {   //### switched from 81 - 100 to have passive cells
	noisetime=0
	gidnoise= i+ ii*N 
    if(pnm.gid_exists(gidnoise)) {
	nc=new NetCon(nil, pnm.pc.gid2obj(gidnoise).synlist.o(3))
	//print gidnoise, pnm.pc.gid2obj(gidnoise).synlist.o(3), nc.precell, nc.postcell, nc
	nc.weight=.00001  //could not compare with the synapses.par file because those are corrected by .001 or (1e-5 x area)
	noiselist.append(nc)
	coweight=nc.weight
	while(noisetime<tStop){
		if (noisetime>0 && noisetime<tStop) { 
		  rng=noiseR.repick()
		  if (rng>basketthreshold) {
			amp=noiseA.repick()
			
			noiselist.o[i].weight = nc.weight * amp
			print noiselist.o[i], noisetime, noiselist.o[i].weight 
			sprint(cmd, "%s.active(1)",noiselist.o[i])
			cvode.event(noisetime,cmd)
			//print i, noisetime
			noiselist.o[i].event(noisetime)

			} //rng
			
			nc.weight=coweight  //reset weight to original value, else it retains it from this event
		  } //if noisetime
		

		noisetime=noisetime+noisestepper
		} //while noisetime

    } //gidexists
	}	//for i


} //ii
//print noiselist.o[81].postcell, noiselist.o[81].weight, noiselist.o[81]
		

} // proc randomPyrThenStim


// =================================================================================================
//
// activeSynapsesRandom()     This is the main function for producing noise, both in SR and CR simulations
//
// =================================================================================================


proc activeSynapsesRandom() {localobj xo, co, noiseR, noiseA
//$1 is tstop, $2 is pyr noise threshold, $3 is bask noise threshold
	synR = new Random( startsw() )

noiseR = new Random( startsw() )
noiseA = new Random( startsw() )
	synR.uniform(5,15) 	
	co   = pnm.nclist
//	if ($1) print "synapse ACTIVATION at ", pnm.pc.id, co.count()
	for ltr(xo, co) {		
		sprint(cmd, "%s", xo.postcell)		
		if(strcmp(cmd, "NULLobject")!=0) {			
			sprint( cmd, "%s.active(1)", xo )			
			cvode.event(synR.repick(), cmd)    //makes NetCon[#] active at random times within 250, 750
   	
		}
	}
//add a noise event to every pyr cell  and every bask cell
noiseR.uniform(0,100)
noiseA.poisson(0.8)
threshold=$2  //75 makes it start out at 100 hz
basketthreshold=$3
noisestepper=.1  //how long the time intervals are in ms between each poisson noise iteration
rng=0
tStop=$1
noiselist=new List()
	noisetime=0
// set up the Netcons for noise synapses on pyramidal cells
for ii=0, nModules-1 {
   for i=0,80 {
	noisetime=0   //will make a delay before noise starts
	gidnoise= i+ ii*N
    if(pnm.gid_exists(gidnoise)) {
		nc=new NetCon(nil, pnm.pc.gid2obj(gidnoise).synlist.o(2))
	nc.weight=.00001 * 53.407075   //cell area of pyr vs. the bask
	noiselist.append(nc)
//###this is a potential race problem--with multiple nodes, will the noiselist number still correspond to gid??
	coweight=nc.weight

	//start=noiseR.repick()
	//start=start/20    // to make scattered delay starts, so they aren't all at same time exactly
	start=0

	while(noisetime<tStop){
		if (noisetime>start && noisetime<tStop) {   //method is redundant, was used for variable input timing
		  rng=noiseR.repick()
		  if (rng>threshold) {
			amp=noiseA.repick()
			noiselist.o[i].weight = nc.weight * amp
			//print noiselist.o[i], noisetime, noiselist.o[i].weight 
			sprint(cmd, "%s.active(1)",noiselist.o[i])
			cvode.event(noisetime,cmd)
			//print i, noisetime
			noiselist.o[i].event(noisetime)
			} //rng			
			nc.weight=coweight  //reset weight to original value, else it retains it from this event
		  } // if noisetime
		noisetime=noisetime+noisestepper
		} //while noisetime
    } //gidexists
	} // for i
  for i=81,100 {   
	noisetime=0
	gidnoise= i+ ii*N
    if(pnm.gid_exists(gidnoise)) {
	nc=new NetCon(nil, pnm.pc.gid2obj(gidnoise).synlist.o(3))
	nc.weight=.00001  //could not compare with the synapses.par file because those are corrected by .001 or (1e-5 x area)
	noiselist.append(nc)
	coweight=nc.weight
	while(noisetime<tStop){
		if (noisetime>0 && noisetime<tStop) { 
		  rng=noiseR.repick()
		  if (rng>basketthreshold) {
			amp=noiseA.repick()			
			noiselist.o[i].weight = nc.weight * amp
			sprint(cmd, "%s.active(1)",noiselist.o[i])
			cvode.event(noisetime,cmd)
			//print i, noisetime
			noiselist.o[i].event(noisetime)
			} //rng			
			nc.weight=coweight  //reset weight to original value, else it retains it from this event
		  } //if noisetime
		noisetime=noisetime+noisestepper
		} //while noisetime
    } //gidexists
  }	//for i
} //ii

		

}  //proc

// =================================================================================================
//
// signalSynapse()
//
// =================================================================================================

proc signalSynapse() {localobj xo, co
//$1 is tstop, $2 is synapse type: 0 for AMPA; 1 for NMDA, $3 is freq in Hz, $4 is the gid getting the sig
	co   = pnm.nclist
tStop=$1
syntype=$2
freq=$3
gidsig=$4
if (freq==0) return
	signaltime=0
// set up the Netcons for signal synapses on pyramidal cells
	if(pnm.gid_exists(gidsig)) {
		if (syntype==0) {
			nc=new NetCon(nil, pnm.pc.gid2obj(gidsig).synlist.o(3))
		}else if (syntype!=0) {
			print "ERROR: do not have any signal other than AMPA now"
			return
		}
	nc.weight=.0001 * 53.407075   //cell area of pyr vs. the bask 
//.0001 *53.4 is still subthreshold .0003 is just right for over threshold
	signallist.prepend(nc)
	while(signaltime<tStop){		
			sprint(cmd, "%s.active(1)",signallist.o[0])
			cvode.event(signaltime,cmd)
			signallist.o[0].event(signaltime)
			signaltime= signaltime + 1000/freq  //to make it ms
		
		} //while noisetime
    } //gidexists

}  //proc



proc noiseStim() {
print "noiseStim is not functional!"
}

// =================================================================================================
//
// addGap(g, gid)  this function is not working, it is done in a different fashion
//
// =================================================================================================

proc addGap() {
	print "addGap does not work! "
	pnm.pc.gid2obj($2).addGapA1($1)
	// all are set to the number listed in createPyrCells now

}

strdef modFileName
// =================================================================================================
//
// addSyn(synParam, gid)
//
// =================================================================================================
proc addSyn() { local i, tao1, tao2, Erev, synLocSec, synLoc, gid, ind, synID, r localobj synParamSet

	synParamSet = $o1
	gid         = $2
	
	for i=0, synParamSet.synSet.count()-1 {
		if (!strcmp(synParamSet.synSet.object(i).postCell, $s3)) {		
			tao1        = synParamSet.synSet.object(i).tao1
			tao2        = synParamSet.synSet.object(i).tao2
			Erev        = synParamSet.synSet.object(i).Erev
			modFileName = synParamSet.synSet.object(i).modFileName
			synLocSec   = synParamSet.synSet.object(i).synLocSec
			synLoc      = synParamSet.synSet.object(i).synLoc
			synID       = synParamSet.synSet.object(i).synID
			r           = synParamSet.synSet.object(i).r
//print synParamSet.synSet.object(i).preCell ," ", synID," ",synParamSet.synSet.object(i).postCell , i
//first time through synID is all -2, then become 0,1,2,3. The precells::i are-- Bask::0 OLM::3 Noise::7 Pyr::9
//
//on the baskets, there are also 0,1,2,3 a bask 1, NMDA pyr 2, and OLM 4, and noise 8
//and on OLM there are 0,1, bask and NMDA pyr
//
			
			// Assign synID
			if (synLocSec==0) {
       			synParamSet.synSet.object(i).synID = pnm.pc.gid2obj(gid).addSynS( tao1, tao2, Erev, modFileName, synLoc)
				}
			
			if (synLocSec==-1) {
       			synParamSet.synSet.object(i).synID = pnm.pc.gid2obj(gid).addSynB( tao1, tao2, Erev, modFileName, synLoc)
				}
			
			if (synLocSec==1) {
       			synParamSet.synSet.object(i).synID = pnm.pc.gid2obj(gid).addSynA1( tao1, tao2, Erev, modFileName, synLoc)
				}
			
			if (synLocSec==2) {
       			synParamSet.synSet.object(i).synID = pnm.pc.gid2obj(gid).addSynA2( tao1, tao2, Erev, modFileName, synLoc)
				}
			
			if (synLocSec==3) {
       			synParamSet.synSet.object(i).synID = pnm.pc.gid2obj(gid).addSynA3( tao1, tao2, Erev, modFileName, synLoc)
				}
			
			if(r>-1) pnm.pc.gid2obj(gid).synlist.object(synParamSet.synSet.object(i).synID).r = r
			
			
			if(synParamSet.synSet.object(i).synID<0) print "ERROR!!!!!!!!*********************"
		}
	}
}


// =================================================================================================
//
// recordVoltages()
//
// =================================================================================================
proc recordVoltages() { local i, ii
	
	for i=0,Npyr-1 {
		if (pnm.gid_exists(i)) pnm.pc.gid2obj(i).recordVoltage()
	}
}




// =================================================================================================
//
// writeVoltages()              1/31/08 I am altering this section
//
// =================================================================================================
proc writeVoltages() { local i, ii  localobj fo, ve, m, mm   //added the localobj
	
	mm= new Matrix(pnm.pc.gid2obj(0).recordT[0].size(), 2)   //I added, not (if.existed)
	ve = new Vector(pnm.pc.gid2obj(0).recordT[0].size())

	mm.setcol(0, pnm.pc.gid2obj(0).recordT[0])
	for i=0,Npyr-1 {
		if (pnm.gid_exists(i)) {			
				ve.add(pnm.pc.gid2obj(i).voltageRecS[1])	
			}
		}     //adds voltages from all pyramidal cells together
	ve.mul(1/(Npyr-1))   //then averages them
	mm.setcol(1, ve)   //then sets them into the matrix that will be saves as the "sum" file
			
	fo = new File()
	if (numarg()==4)	{sprint(cmd, "data/sum_b%d_p%d_g%d_f%d.dat",$1,$2,$3,$4) //3.1lf for i nA
	} else {sprint(cmd, "data/sum.dat")}
	
	fo.wopen(cmd)
	mm.fprint(0, fo, "%6.3lf ")
	fo.close()
}

// =================================================================================================
//
// writesingleVoltage()              12/18/08 for when there is just a single pyr cell conn
//
// =================================================================================================
proc writesingleVoltage() { localobj fo,  mm   

//$1 is the basketthr $2 is pyrthresh $3 is sigfreq
//I am assuming the cell we are recording is "GID= 1"

fo= new File()
	if (pnm.gid_exists(1)){
		mm= new Matrix(pnm.pc.gid2obj(1).recordT[1].size(), 2)   
		mm.setcol(0, pnm.pc.gid2obj(1).recordT[1])
		mm.setcol(1, pnm.pc.gid2obj(1).voltageRecS[1])
     	
		sprint(cmd, "data/1_b%d_p%d_single_f%d.dat",$1,$2,$3) 
		fo.wopen(cmd)
		mm.fprint(0, fo, "%6.3lf ")
		fo.close()

	}
	
   	if (pnm.gid_exists(81)) {
		mm= new Matrix(pnm.pc.gid2obj(81).recordT.size(), 2)   
		mm.setcol(0, pnm.pc.gid2obj(81).recordT)
		mm.setcol(1, pnm.pc.gid2obj(81).voltageRecS)
     	
		sprint(cmd, "data/81_b%d_p%d_single_f%d.dat",$1,$2,$3) 
		fo.wopen(cmd)
		mm.fprint(0, fo, "%6.3lf ")
		fo.close()
	}
}


proc recordnoise() { 
	pnm.pc.gid2obj(1).recordnoise()
}

// =================================================================================================
//
// writenoise()              12/18/08 to determine noise input variance
//
// =================================================================================================
proc writenoise() { localobj fo,  mm   

//$1 is the pyrthresh 
//I am assuming the cell we are recording is "GID= 1"

	fo= new File()
	if (pnm.gid_exists(1)){
		mm= new Matrix(pnm.pc.gid2obj(1).recordT[1].size(), 2)   
		mm.setcol(0, pnm.pc.gid2obj(1).recordT[1])
		mm.setcol(1, pnm.pc.gid2obj(1).iRec)
     	
		sprint(cmd, "data/noise_b%d.dat",$1) 
		fo.wopen(cmd)
		mm.fprint(0, fo, "%6.3lf ")
		fo.close()

	}
}


// =================================================================================================
//
// setScatteredVoltages(low, high)
//
// =================================================================================================
proc setScatteredVoltages() { local low, high

  low  = $1
  high = $2

  i   = 0
  gid = 0
  
  for pyrIter(&i,&gid)  if (pnm.gid_exists(gid)) pnm.pc.gid2obj(gid).setScatteredVoltages(low, high)  
//for pyrIter(&i,&gid)  if (pnm.gid_exists(gid)) pnm.pc.gid2obj(gid).setScatteredVoltages(-56, -55)
  for baskIter(&i,&gid) if (pnm.gid_exists(gid)) pnm.pc.gid2obj(gid).setScatteredVoltages(low, high)
  for olmIter(&i,&gid)  if (pnm.gid_exists(gid)) pnm.pc.gid2obj(gid).setScatteredVoltages(low, high)


  finitialize() 
 
  finitialize() 
//if (debugFlag) print "not stuck there"
}





// =================================================================================================
//
// makePyrgCaScattered( flag ) flag = {0,1}
//
// =================================================================================================
proc makePyrgCaScattered() { local i, low, high, gid localobj rand

	gCaScattered = $1	
	gid = 0
	i   = 0
	rand = new Random(startsw())
	if (gCaScattered==0) {	
		low  = 10
		high = 10		
	}else{	
		low  =  9
		high = 11	
	}	
	rand.uniform(low, high)
    for pyrIter(&i, &gid) if (pnm.gid_exists(i)) pnm.pc.gid2obj(i).den.gca_icapr = rand.repick()
}



strdef lineStr

// =================================================================================================
//
// connectNetwork()
//
// =================================================================================================
proc connectNetwork() { local sgid, i, ii, thr, Erev, g, cellArea, num, cnt, countSyn, preGID, postGID, noiseGID, synID, globalID, ind, gmax, gmaxUnits, delay, Npre, gapcon localobj fo
		
	preGID   = 0
	postGID  = 0
	noiseGID = 0
	globalID = 0
	cellArea = 0
	cnt=0

    Npyr  = synCellSet.cellSet.object(0).N + 1  //I think this is not supposed to have +1 for 80 cells
    Nbask = synCellSet.cellSet.object(1).N
    Nolm  = synCellSet.cellSet.object(2).N
   par_gaps=new List()    //makes a gap junction list

   for ii=0, nModules-1 {
	fo = new File()
	fo.ropen("parameters/conn.dat")
	countSyn = 0
	countGap = 0

	while( !fo.eof() ) {
	
		fo.gets( lineStr )
		
		cnt=cnt+1
		num = sscanf( lineStr, "%d, %d, %d", &preGID, &postGID, &globalID)
		preGID  = preGID  + ii * N
		postGID = postGID + ii * N    // has GABA (globalID 0) from each bask to pyr once
							// then has random NMDA pyr to bask (GlobalID 2, which is NOT synID)
							// these are the corresponding 0 and 2 rows of the synapses.par file
							// I have now made GlobalID 7 the noise->pyr AMPA,  and 8 noise-> bask
							// 9 signal->pyr 10 pyr-pyr gaps, 11 bask-bask gaps
		
		if(num!=3) print "ERROR!!!!!!!!!!!!!!"		
		//gAMPAdend * taPyr * 1e-5

		// Do only if source or target are on the node
		if(pnm.gid_exists(preGID) || pnm.gid_exists(postGID)) {

			gmax  = synParamSet.synSet.object(globalID).gmax
			delay = synParamSet.synSet.object(globalID).delay
			synID = synParamSet.synSet.object(globalID).synID
					// these get set in AddSyn
			gmaxUnits = synParamSet.synSet.object(globalID).gmaxUnits
			Erev  = synParamSet.synSet.object(globalID).Erev
			if ((manParamFlag==1) && (globalID==0) ) {
					print "manual param flag ", preGID, postGID
					//gmax=	gmax_basktopyr
				}


			if ( pnm.gid_exists(postGID)) cellArea = pnm.pc.gid2cell(postGID).getTotalArea()
			
			if(gmaxUnits==1) {
			
				// Convert gmax from mS/cm2 to uS
				g = cellArea * gmax * 1e-5   //cellArea is 100 in bask, so g=gmax*.001
								
			}else{
			
				// Convert gmax from nS to uS
				g = gmax * 1e-3
			}
			
			// Normalize by the presynaptic number of cells
			
			Npre = synParamSet.synSet.object(globalID).Npre
			if (Npre==-1) {   //not sure what this does, it affects those syns in synapses.par with -1 Npre

				pN = -1
				pNN = 0				
				for i=0,  synCellSet.cellSet.count()-1 {			
					if(strcmp(synCellSet.cellSet.object(i).name, synParamSet.synSet.object(globalID).preCell)==0) {					
						pN = synCellSet.cellSet.object(i).N
						pNN = pNN + 1
					}			
				}				
				if(pNN!=1 || pN==-1) print "ERROR pNN"
				Npre = pN				
			}			
			g = g / Npre
			
				// the gmax is divided by 1000, then by the presyn cells (20 in case of bask-pyr GABA)
			if (globalID<10) {
				if(globalID==0){    //this is for the basket-pyr AMPA synapses, which I can split into 2 (but did not in this paper)
			         {   //if(cnt%2)     ###this is to make only every other one go if placed before the bracket
				    ind   = pnm.nc_append(preGID, postGID, synID, g, delay) 
					
			        }
                        } else {ind   = pnm.nc_append(preGID, postGID, synID, g, delay) 
					
					}
			
			     } else {	

			// globalID 10 and 11 are gaps.  we've already read in the gmax of the gap, rest is garbage

				if (globalID == 10) { //for baskbask gaps
					
					par_gap_create(preGID, 0 ,postGID, 0, gmax )   //the 2nd and 4th are locations
					countGap = countGap+1					// 0 means soma, 1 means A1 for now
				}  //globalID 10
				if (globalID == 11) { //for pyrpyr gaps
					if (manParamFlag==1) {
						print "manual param flag ", preGID, postGID
						//gmax=	gmax_pyrgap
	
					}  //manParamFlag

					par_gap_create(preGID, 0 ,postGID, 0, gmax )   //the 2nd and 4th are locations
					countGap = countGap+1					// 0 means soma, 1 means A1 for now
				}  //globalID 11
			} //else
			
			countSyn = countSyn + 1
			thr=pnm.pc.threshold(preGID,-20)  //can add (preGID, -10) to change it. At 0, lost many conns
		}//if pre/post GIDs exist
				 
	} //while f.!eof
    
	fo.close()

   } // ii


countSyn=countSyn-countGap //to correct for increments of Syn that were actually Gap above
totalGaps=totalGaps+countGap
ggap=.00005
	pnm.pc.setup_transfer()
	printf("%d synapses and %d gaps were set\n\n", countSyn, totalGaps)
	
}



// =================================================================================================
//
// connectNetworkInterModules()     THIS SECTION NOT UTILIZED IN THIS MODEL
//
// =================================================================================================
proc connectNetworkInterModules() { local sgid, i, ii, iii, Erev, g, cellArea, num, countSyn, preGID, postGID, synID, globalID, ind, gmax, gmaxUnits, delay, Npre localobj fo
	//countSyn = 0	//for some reason it's passing the N as the countSyn now 5/8 WS
	preGID   = 0
	postGID  = 0
	globalID = 0
	cellArea = 0
	
	for ii=0, nModules-1 {
	 for iii=0,nModules-1 if (ii!=iii) {

	  fo = new File()
	  fo.ropen("parameters/conn.dat")
	  countSyn = 0
	
	
	   while( !fo.eof() ) {
	
		fo.gets( lineStr )
		
		num = sscanf( lineStr, "%d, %d, %d", &preGID, &postGID, &globalID)
		
		if(strcmp(synParamSet.synSet.object(globalID).preCell, "OLM")==0 && strcmp(synParamSet.synSet.object(globalID).postCell, "Pyr")==0 && preGID!=postGID) {
		
		 preGID  = preGID  + ii * N
		 postGID = postGID + iii * N
		 if(num!=3) print "ERROR!!!!!!!!!!!!!!"		
		 //gAMPAdend * taPyr * 1e-5

		 // Do only if source or target are on the node
		  if(pnm.gid_exists(preGID) || pnm.gid_exists(postGID)) {

			gmax  = synParamSet.synSet.object(globalID).gmax
			delay = synParamSet.synSet.object(globalID).delay
			synID = synParamSet.synSet.object(globalID).synID
			gmaxUnits = synParamSet.synSet.object(globalID).gmaxUnits
			Erev  = synParamSet.synSet.object(globalID).Erev
			
			if ( pnm.gid_exists(postGID)) cellArea = pnm.pc.gid2cell(postGID).getTotalArea()
			
			if(gmaxUnits==1) {
			
				// Convert gmax from mS/cm2 to uS
				g = cellArea * gmax * 1e-5
				
			}else{
			
				// Convert gmax from nS to uS
				g = gmax * 1e-3
		  }
			
			// Normalize by the presynaptic number of cells
			
			Npre = synParamSet.synSet.object(globalID).Npre
			
			if (Npre==-1) {

				pN = -1
				pNN = 0
				
				for i=0,  synCellSet.cellSet.count()-1 {
		
					if(strcmp(synCellSet.cellSet.object(i).name, synParamSet.synSet.object(globalID).preCell)==0) {
						
						pN = synCellSet.cellSet.object(i).N
						pNN = pNN + 1
					}			
				}
				
				if(pNN!=1 || pN==-1) print "ERROR pNN"
				Npre = pN
				
			}
			
			g = g / Npre
			
			ind   = pnm.nc_append(preGID, postGID, synID, g, delay)
			
			countSyn = countSyn + 1
		}
	}
    }
	fo.close()
	} // ii

	}// iii

	pnm.pc.setup_transfer()
	printf("%d Intermodule Connections were set\n\n", countSyn)

}
// =================================================================================================
//
// par_gap_create()
//  taken from NEURON implementation of Traub's 2005 thalamocortical model
// =================================================================================================

gap_src_gid = 2
func par_gap_create() {
	gap_src_gid += 2

	if (pnm.gid_exists($1)) { 
		par_gap_create1($1, $2, gap_src_gid + 1, gap_src_gid, $5)
	}
	if (pnm.gid_exists($3)) {
		par_gap_create1($3, $4, gap_src_gid, gap_src_gid + 1, $5)
	}
	return totalGaps + 1
}
proc par_gap_create1() {localobj c, g, gaploc
	c = pnm.pc.gid2obj($1)
	if ($2==0) {gaploc=c.getlocoS()}
	if ($2==1) {gaploc=c.getlocoA1()}
	gaploc.sec {      
			g = new gGapPar(.5)
			par_gaps.append(g)
			//print pnm.pc.gid2obj($1), gaploc, $3, $4, par_gaps.count()
			pnm.pc.target_var(&g.vgap, $3)
			pnm.pc.source_var(&v(.5), $4)
			g.g = $5		
	}
}



endtemplate TGnet