Network recruitment to coherent oscillations in a hippocampal model (Stacey et al. 2011)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:135903
"... Here we demonstrate, via a detailed computational model, a mechanism whereby physiological noise and coupling initiate oscillations and then recruit neighboring tissue, in a manner well described by a combination of Stochastic Resonance and Coherence Resonance. We develop a novel statistical method to quantify recruitment using several measures of network synchrony. This measurement demonstrates that oscillations spread via preexisting network connections such as interneuronal connections, recurrent synapses, and gap junctions, provided that neighboring cells also receive sufficient inputs in the form of random synaptic noise. ..."
Reference:
1 . Stacey WC, Krieger A, Litt B (2011) Network recruitment to coherent oscillations in a hippocampal computer model. J Neurophysiol 105:1464-81 [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
//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


//this template is very similar to that in the 2009 paper: ModelDB accession 135902
//it has the following changes:
//1. addition of 20 more pyramidal cells "Neighbors", taking the gid slots 81-100.  Original 0-80 are "Drivers"
//2. basket cells now take gid slots 101-120
//3. basket cells contact all 100 pyr cells (GABA), and each pyr cell contacts 2 or 3 basket cells (AMPA/NMDA)
//4. potential recurrent synapse or gap junction conn between 20 of the Driver cells and the 20 Neighbors
//5. Only inputs are noise to the pyramidal cells, implemented in the procedure activeSynapsesRandom()
//    that function starts with low noise for 200 ms, then adds noise to the Drivers for 1000 ms, then at 600
//    ms adds noise to the Neighbors for 1000 ms (thus 600 ms crossover).  Total duration 1600 ms
//6. In the writeVoltages file, now saves three average voltage files: for all pyr cells ("sum")
//     all Driver cells "sumdrive" and all Neighbor cells "sumpas"


begintemplate Recruitnet

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("../parameters/synapses.tem")}
{load_file("../parameters/recruitcells.tem")}



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

	oneCell = 0   //I don't know 
	nModules  = 1    // this was 4 in the Tort paper
	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()
	//print N, Npyr, Nbask
	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)
	}
}


// =================================================================================================
//
// manualParams()   $1 is the gmax for the bask-pyr GABA, $2 is gmax for pyr-pyr gap junction
//
// =================================================================================================
proc manualParams() {local i
	manParamFlag=1
	gmax_basktopyr=$1
	gmax_pyrgap=$2
	print manParamFlag, "entered manualparams, not yet implemented"

}

// =================================================================================================
//
// 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 )
		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"
		// Random vurrent inj 0-500 ms
		// c.iappS.set_random(gid)		
		//c.recordVoltage()
    }
}

// =================================================================================================
//
// createNoiseStim()   No longer used
//
// =================================================================================================
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


		//addSyn(synParamSet, gid, "OLM")
		
		//print "creating noise cells"
		
    }
}



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.  THIS WAS NOT USED FOR THIS PAPER
// =================================================================================================


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 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, as basket noise was treated in the 2009 paper
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()

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

   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 to add noise to the network
//    It has been altered from the 2009 model significantly
// =================================================================================================

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

	noiseR = new Random( startsw() )
	noiseA = new Random( startsw() )
	
	synR.uniform(5,15)
 	co   = pnm.nclist
	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
   
			//###!  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)
	threshold=$2  //75 makes it start out at 100 hz
	basketthreshold=$3
	driverthreshold=$4
	noisestepper=.1  //how long the time intervals are in ms between each poisson noise iteration
	delayD=200 //how long to wait before the Driver noise turns on
	lengthD=1000 //how long the Driver noise will last
	delayP=600 //how long for the Neighbor cells to wait before noise starts, now connected to "start" below
	quiescentthr=97 //baseline low noise level before it starts up higher at delayD

	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,Npyr-1 {
	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)
		coweight=nc.weight

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

     if(gidnoise<81){    //for the Driver Cells
	  while(noisetime<(lengthD+delayD)){
		if (noisetime<delayD) {   //what to do before the increased driver noise starts
		  rng=noiseR.repick()
		  if (rng>quiescentthr) {
			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  //reset weight to original value, else it retains it from this event
		  }  else  {   //after the delay  
		  rng=noiseR.repick()
		  if (rng>driverthreshold) {
			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  //reset weight to original value, else it retains it from this event
		  } // if noisetime
		noisetime=noisetime+noisestepper
	  } //while noisetime

     } else {   //the 20 Neighbor cells with different threshold

	while(noisetime<tStop){
		if (noisetime<start) {
		  rng=noiseR.repick()
		  if (rng>quiescentthr) {
			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  //reset weight to original value, else it retains it from this event
		  } 	else {   //after the off period
		//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 
			sprint(cmd, "%s.active(1)",noiselist.o[i])
			cvode.event(noisetime,cmd)
			noiselist.o[i].event(noisetime)
			} //rng
			nc.weight=coweight  //reset weight to original value, else it retains it from this event
		  } // else
		noisetime=noisetime+noisestepper
		} //while noisetime
     } //else
    } //gidexists
   } // for i

   for i=Npyr,N-1 {   //### Basket cells

	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)
			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 ii
}  //proc

// =================================================================================================
//
// signalSynapse()   Not used in this paper
//
// =================================================================================================

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 nonfunctional"
}

// =================================================================================================
//
// 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
//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()     This now saves several files.  "sum" is the sum of all 100 pyr cells         
//				"sumdrive" is the sum of just the Driver pyr cells
//				"sumpas" is the sum of just the Neighbor pyr cells
// =================================================================================================
proc writeVoltages() { local i, ii  localobj fo, ve, vedrive, vepassive, m, mm, mmdrive, mmpassive   //added the localobj

	if (pnm.gid_exists(0)) {
		mm= new Matrix(pnm.pc.gid2obj(0).recordT[0].size(), 2)   
		mmdrive= new Matrix(pnm.pc.gid2obj(0).recordT[0].size(), 2)   
		mmpassive= new Matrix(pnm.pc.gid2obj(0).recordT[0].size(), 2)
		ve = new Vector(pnm.pc.gid2obj(0).recordT[0].size())
		vedrive = new Vector(pnm.pc.gid2obj(0).recordT[0].size())
		vepassive = new Vector(pnm.pc.gid2obj(0).recordT[0].size())  
		mm.setcol(0, pnm.pc.gid2obj(0).recordT[0])
		mmdrive.setcol(0, pnm.pc.gid2obj(0).recordT[0])
		mmpassive.setcol(0, pnm.pc.gid2obj(0).recordT[0]) 
	}

//first make the sum of all pyr cells
	
	for i=0,Npyr-1 {
		if (pnm.gid_exists(i)) {				
				ve.add(pnm.pc.gid2obj(i).voltageRecS[1])					
						}
		}
	ve.mul(1/(Npyr-1))
	mm.setcol(1, ve)		
	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()

//then make the sum of the Driver "drive" cells, and the Neighbor "pas" cells

  if (1) {    //can just switch this to zero if I don't want extra files.	
	for i=0,80 {
		if (pnm.gid_exists(i)) {				
				vedrive.add(pnm.pc.gid2obj(i).voltageRecS[1])	
			}
		}
	vedrive.mul(1/81)
	mmdrive.setcol(1, vedrive)			
	fo = new File()
	if (numarg()==4)	{sprint(cmd, "data/sumdrive_b%d_p%d_g%d_f%d.dat",$1,$2,$3,$4) //3.1lf for i nA
	} else {sprint(cmd, "data/sumdrive.dat")}	
	fo.wopen(cmd)	
	mmdrive.fprint(0, fo, "%6.3lf ")	
	fo.close()
	for i=81,100 {
		if (pnm.gid_exists(i)) {				
				vepassive.add(pnm.pc.gid2obj(i).voltageRecS[1])	
			}
		}
	vepassive.mul(1/20)
	mmpassive.setcol(1, vepassive)			
	fo = new File()
	if (numarg()==4)	{sprint(cmd, "data/sumpas_b%d_p%d_g%d_f%d.dat",$1,$2,$3,$4) //3.1lf for i nA
	} else {sprint(cmd, "data/sumpas.dat")}	
	fo.wopen(cmd)	
	mmpassive.fprint(0, fo, "%6.3lf ")	
	fo.close()
   } //if (1)

}

// =================================================================================================
//
// writesingleVoltage()              12/18/08 for when there is just a single pyr cell conn
//						Not used in this paper
// =================================================================================================
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
//					not used in this paper
// =================================================================================================
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/recruitconn.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
			
			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
  pnm.pc.setup_transfer()
  printf("%d synapses and %d gaps were set\n\n", countSyn, totalGaps)
	
}



// =================================================================================================
//
// connectNetworkInterModules()
//    THIS IS NOT USED
// =================================================================================================
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/recruitconn.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 Recruitnet