Hippocampal CA1 NN with spontaneous theta, gamma: full scale & network clamp (Bezaire et al 2016)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:187604
This model is a full-scale, biologically constrained rodent hippocampal CA1 network model that includes 9 cells types (pyramidal cells and 8 interneurons) with realistic proportions of each and realistic connectivity between the cells. In addition, the model receives realistic numbers of afferents from artificial cells representing hippocampal CA3 and entorhinal cortical layer III. The model is fully scaleable and parallelized so that it can be run at small scale on a personal computer or large scale on a supercomputer. The model network exhibits spontaneous theta and gamma rhythms without any rhythmic input. The model network can be perturbed in a variety of ways to better study the mechanisms of CA1 network dynamics. Also see online code at http://bitbucket.org/mbezaire/ca1 and further information at http://mariannebezaire.com/models/ca1
References:
1 . Bezaire MJ, Raikov I, Burk K, Vyas D, Soltesz I (2016) Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent CA1 circuit. Elife [PubMed]
2 . Bezaire M, Raikov I, Burk K, Armstrong C, Soltesz I (2016) SimTracker tool and code template to design, manage and analyze neural network model simulations in parallel NEURON bioRxiv
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 cell; Hippocampus CA1 interneuron oriens alveus cell; Hippocampus CA1 basket cell; Hippocampus CA1 stratum radiatum interneuron; Hippocampus CA1 bistratified cell; Hippocampus CA1 axo-axonic cell; Hippocampus CA1 PV+ fast-firing interneuron;
Channel(s): I Na,t; I K; I K,leak; I h; I K,Ca; I Calcium;
Gap Junctions:
Receptor(s): GabaA; GabaB; Glutamate; Gaba;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; NEURON (web link to model);
Model Concept(s): Oscillations; Methods; Connectivity matrix; Laminar Connectivity; Gamma oscillations;
Implementer(s): Bezaire, Marianne [mariannejcase at gmail.com]; Raikov, Ivan [ivan.g.raikov at gmail.com];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal cell; Hippocampus CA1 interneuron oriens alveus cell; GabaA; GabaB; Glutamate; Gaba; I Na,t; I K; I K,leak; I h; I K,Ca; I Calcium; Gaba; Glutamate;
/***********************************************************************************************
I.  LOAD LIBRARIES
***********************************************************************************************/
{load_file("nrngui.hoc")}				// Standard definitions - NEURON library file

{load_file("netparmpi.hoc")}			// Contains the template that defines the properties of
										//  the ParallelNetManager class, which is used to set
										//   up a network that runs on parallel processors
										
{load_file("./setupfiles/ranstream.hoc")}	// Contains the template that defines a RandomStream
											//  class used to produce random numbers
											// 	for the cell noise (what type of noise?)
											
{load_file("./setupfiles/CellCategoryInfo.hoc")}	// Contains the template that defines a 
													//  CellCategoryInfo class used to store
													// 	celltype-specific parameters

{load_file("./setupfiles/defaultvar.hoc")}	// Contains the proc definition for default_var proc
													//  that can't be changed at command line

{load_file("./setupfiles/parameters.hoc")}	// Loads in operational and model parameters that can
											//  be changed at command line											
{load_file("./setupfiles/set_other_parameters.hoc")}// Loads in operational and model parameters

/**********************/
default_var("BasedOnRunName","ReserveTrestles_01")		// Name of simulation run that the netclamp is based on
default_var("CellOI","axoaxoniccell")		// type of presynaptic cell
default_var("CellOItype",0)		// type of presynaptic cell
default_var("gid",11)		// gid of presynaptic cell

strdef cmdstr, cmd, RunName

sprint(RunName,"%s_%d", CellOI, gid)
/***********************************************************************************************
II. SET MODEL SIZE, CELL DEFINITIONS
***********************************************************************************************/
celsius=34

{load_file("./setupfiles/load_cell_category_info.hoc")}	// Reads the 'cells2include.hoc' file and
														//  loads the information into one
														//  'CellCategoryInfo' object for each cell
														//  type (bio or art cells?). Info includes
														//  number of cells, gid ranges, type name 

strdef tempFileStr
sprint(tempFileStr,"./cells/class_%s.hoc",cellType[CellOItype].technicalType)	// Concatenate the
																				//  path and file
print "loading ", tempFileStr																
load_file(tempFileStr)			// Load the file with the template that defines the classs

sprint(tempFileStr,"./cells/class_%s.hoc","ppspont")	// Concatenate the
print "loading ", 		tempFileStr																
load_file(tempFileStr)			// Load the file with the template that defines the class

totalCells = 1						// totalCells counts the number of 'real' cells
									
// Create numStims number of NetStims
strdef fn
objref f2
f2 = new File()
sprint(fn, "./netclamp/%s/%s_%d_NetStimConns.dat", BasedOnRunName, CellOI,  gid)
f2.ropen(fn)
numStims = f2.scanvar

ncell = 1 + numStims						// ncell counts all 'real' and 'artificial' cells
							
objref TMPcellType[numCellTypes+1]
								
for i=0,numCellTypes-1 {
	if (i==CellOItype) {
		cellType[CellOItype].numCells = 1	
	} else {
		cellType[i].numCells = 0
	}
	if (i>0) {
		cellType[i].updateGidRange(cellType[i-1].cellEndGid+1)	// Update the gid range for each
	} else {
		cellType[i].updateGidRange(0)	// Update the gid range for each
	}
	TMPcellType[i] = cellType[i]
}

TMPcellType[numCellTypes] = new CellCategoryInfo(numCellTypes)	// Make one object for each cell type to store cell type info
TMPcellType[numCellTypes].setCellTypeParams("PatternStim", "ppspont", 1, numStims, 0, 1)	// Set parameters
TMPcellType[numCellTypes].numCons = new Vector(numCellTypes,0)

numCellTypes = numCellTypes+1

objref cellType[numCellTypes]

for i=0,numCellTypes-1{
	cellType[i] = TMPcellType[i]
	print cellType[i].cellType_string, " ", cellType[i].numCells, "(|", cellType[i].cellStartGid, "| to |", cellType[i].cellEndGid, "|)"
}


objref TMPcellType
/***********************************************************************************************
III.SET UP PARALLEL CAPABILITY AND WRITE OUT RUN RECEIPT
***********************************************************************************************/
objref pnm, pc, nc, nil
proc parallelizer() {
	pnm = new ParallelNetManager(ncell)	// Set up a parallel net manager for all the cells
	pc = pnm.pc
	pnm.round_robin()					// Incorporate all processors - cells 0 through ncell-1
										//	are distributed throughout the hosts
										//	(cell 0 goes to host 0, cell 1 to host 1, etc)
}
parallelizer()


iterator pcitr() {local i2, startgid	// Create iterator for use as a standard 'for' loop
										//  throughout given # cells usage:
										//  for pcitr(&i1, &i2, &gid, it_start, it_end) {do stuff}
										//  it_start and it_end let you define range over
										//  which to iterate
										//  i1 is the index of the cell on the cell list for that host
										//  i2 is the index of that cell for that cell type on that host
	numcycles = int($4/pc.nhost)
	extra = $4%pc.nhost
	addcycle=0
	if (extra>pc.id) {addcycle=1}
	startgid=(numcycles+addcycle)*pc.nhost+pc.id
	i1 = numcycles+addcycle // the index into the cell # on this host.
	i2 = 0 // the index of the cell in that cell type's list on that host
	if (startgid<=$5) {
		for (i3=startgid; i3 <= $5; i3 += pc.nhost) {	// Just iterate through the cells on
														//  this host(this simple statement
														//  iterates through all the cells on
														//  this host and only these cells because 
														//  the roundrobin call made earlier dealt
														//  the cells among the processors in an
														//  orderly manner (like a deck of cards)
				$&1 = i1
				$&2 = i2
				$&3 = i3
				iterator_statement
				i1 += 1
				i2 += 1
		}
	}
}


// Create one Real Cell

objref cells, ransynlist, ranstimlist
cells = new List()						
ransynlist = new List()
ranstimlist = new List()

func xpos_algorithm() { // Arguments: 1-gid, 2-numCells, 3-startGid, 4-binNum, 5- binNum, 6-binSize; Return: x position of cell
	return 0
}
func ypos_algorithm() { // Arguments: gid, numCells, startGid, binNum, binNum, binSize; Return: y position of cell
	return 0
}
func zpos_algorithm() { // Arguments: 1-gid, 2-numCells, 3-startGid, 4-binNum, 5-binSize, cell layer Zo; Return: z position of cell
	return 0
}
													
{load_file("./setupfiles/create_cells_pos.hoc")}	// Creates each cell on its assigned host
													//  and sets its position using the algorithm
													//  defined above



objref cell
for r=cellType[numCellTypes-1].cellStartGid, cellType[numCellTypes-1].cellEndGid {
	cell = pc.gid2cell(r)
	cell.interval=1e9
	cell.start=-1
}

// Add a recording to its soma potential and somatic ion channel currents, as well as a few other locations

strdef nickname, secstr
nickname="soma"
myi_flag=1

strdef cmdstr
strdef mname, tmpstr
objref strobj, mt, cell
objref mechstring[9], mechlength

mechlength = new Vector(1)
strobj = new StringFunctions()

objref cell, mt
cell = pc.gid2cell(cellType[CellOItype].cellStartGid)
secstr="soma"

sprint(cmdstr,"objref mytrace%s", nickname)
{execute1(cmdstr)}

{sprint(cmdstr,"mytrace%s = new Vector(tstop/dt)", nickname)}
{execute1(cmdstr)}

//{sprint(cmdstr,"mytrace%s.record(&cell.%s.v(0.5))", nickname, secstr)}
{sprint(cmdstr,"mytrace%s.record(&%s.v(0.5))", nickname, cell.myroot)}
{execute1(cmdstr)	}	

/* Get mechanisms from soma section "myroot" recorded */

if (myi_flag==1) {
	mt = new MechanismType(0)
	objref mechstring[mt.count()]
	k = 0
	for i=0, mt.count()-1 {
		mt.select( i )
		mt.selected(mname)
		if( ismembrane(mname)) {
			if (strcmp(mname,"capacitance")!=0 && strcmp(mname,"morphology")!=0 && strcmp(mname,"iconc_Ca")!=0  && strcmp(mname,"iconcCa")!=0 && strcmp(mname,"ccanl")!=0 && strcmp(mname,"cad")!=0 && strcmp(mname,"pas")!=0 && strcmp(mname,"vmax")!=0 ) { //
				if (strobj.substr(mname,"_ion")==-1) {
					//printf("myi_%s \n", mname) 
					sprint(tmpstr, "myi_%s", mname) // "cell.soma.%s(0.5)", tmpstr
					mechstring[k] = new String()
					mechstring[k].s = tmpstr
					if (strcmp("myi_ch_Nav", mechstring[k].s)==0) {
					print nickname, ": ", tmpstr}
					{k = k+1}
				} 
			} 
		}
	}
	{mechlength.x[0] = k}
		if (strcmp("som", nickname)==0) {
	print "num mechs: ", k}

	for rt = 0, mechlength.x[0]-1 {
		sprint(cmdstr, "{objref %s_%s}", nickname, mechstring[rt].s)
		execute1(cmdstr)
		sprint(cmdstr, "{%s_%s = new Vector(%g)}", nickname, mechstring[rt].s, (tstop-tstart)/dt)
		execute1(cmdstr)
		sprint(cmdstr, "{%s_%s.record(&%s.%s(0.5))}", nickname, mechstring[rt].s, cell.myroot, mechstring[rt].s)
		execute1(cmdstr)
		if (strcmp("myi_ch_Nav", mechstring[rt].s)==0) {
		print "record: ", cmdstr}
	}
}



// Connect NetStims to all CellOI input synapses
{load_file("./setupfiles/nc_append_functions.hoc")}	// Defines nc_append and nc_appendo, which 
{load_file("./setupfiles/load_cell_conns.hoc")}	// Load in the cell connectivity info

pregid = 1

numConnTypes = f2.scanvar
postcellType = CellOItype

for r=0, numConnTypes-1 {
	precellType = f2.scanvar
	numSynTypes = f2.scanvar
	synWeight = cellType[precellType].wgtConns.x[postcellType]	
		
	for SynNumber = 0, numSynTypes-1 {
		print "pregid = ", pregid, ", precellType = ", precellType, ", synapse = ", SynNumber, ", weight = ", synWeight
		nc_append(pregid, cellType[postcellType].cellStartGid, precellType, SynNumber, synWeight + (pregid+1)*1000, 0.5)	// Make the connection
		pregid += 1
		cellType[postcellType].numCons.x[precellType] += 1
	}
}

f2.close()

// Read in firing times of all net stims
strdef fn
objref f2
objref pattern_, tvec_, idvec_
f2 = new File()
sprint(fn, "./netclamp/%s/%s_%d_NetStimTimes.dat", BasedOnRunName, CellOI,  gid)
f2.ropen(fn)
numspikes = f2.scanvar

pattern_ = new PatternStim()

tvec_ = new Vector(numspikes)
idvec_ = new Vector(numspikes)

for i=0, numspikes-1 {
	tvec_.x[i] = f2.scanvar // spike time in ms
	idvec_.x[i] = f2.scanvar // gid of NetStim to make fire
}
f2.close()

pattern_.fake_output = 1
pattern_.play(tvec_, idvec_)

// Run simulation & write out spike raster and recording results

proc init() { local dtsav, temp, secsav, secondordersav	// Redefine the proc that initializes the
														//  simulation (why?)

	dtsav = dt						// Save the desired dt value to reset it after temporarily
									//  changing dt to run a quick presimulation to allow the
									//  network to reach steady state before we start 'recording'
									
	secondordersav = secondorder	// Save the desired secondorder value to reset it after
									//  temporarily changing secondorder to run a quick presimulation
									//  to allow the network to reach steady state before we start
									//  'recording' its results

	finitialize(v_init)	// Call finitialize from within the custom init proc, just as the default
						//  init proc does. Note that finitialize will call the INITIAL block for
						//  all mechanisms and point processes inserted in the sections and set the
						//	initial voltage to v_init for all sections

	t = -200			// Set the start time for (pre) simulation; -500 ms to allow the network to
						// reach steady state before t = 0, when the real simulation begins
						
	dt= 10				// Set dt to a large value so that the presimulation runs quickly
	
	secondorder = 0		// Set secondorder to 0 to set the solver to the default fully implicit backward
						//  euler for numerical integration (see NEURON ref)
		
	temp= cvode.active()			// Check whether cvode, a type of solver (?) is on
	if (temp!=0) {cvode.active(0)}	// If cvode is on, turn it off while running the presimulation
															
	while(t<-100) { fadvance() if (PrintTerminal>1) {print t}}	// Run the presimulation from t = -500
															//  to t = -100 (why not 0?) to let the
															//  network and all its components reach
															//  steady state. Integrate all section
															//  equations over the interval dt,
															//  increment t by dt, and repeat until
															//  t at -100
															
	if (temp!=0) {cvode.active(1)}	// If cvode was on and then turned off, turn it back on now
	
	t = tstart 						// Set t to the start time of the simulation
	
	dt = dtsav						// Reset dt to the specified value for the simulation
	
	secondorder = secondordersav	// Reset secondorder to the specified value for the simulation
	
	if (cvode.active()){
		cvode.re_init()				// If cvode is active, initialize the integrator
	} else {
		fcurrent()					// If cvode is not active, make all assigned variables
									//	 (currents, conductances, etc) consistent with the
									//   values of the states
	}
}

simnum=0

strdef cmdstr
sprint(cmdstr,"mkdir ./netclamp/%s/%s", BasedOnRunName,RunName)
print cmdstr
system(cmdstr)

proc spikeout() {local i, rank  localobj f				// Write out a spike raster (cell, spike time)
	pc.barrier()									// Wait for all ranks to get to this point
	sprint(cmd,"./netclamp/%s/%s/spikeraster_%03.0f.dat", BasedOnRunName, RunName, simnum)
	f = new File(cmd)
	if (pc.id == 0) { 								// Write header to file 1 time only
		f.wopen()
		f.close()
	}
	
	for rank = 0, pc.nhost-1 {				// For each processor, allow processor to append to file the spike times of its cells
		if (rank == pc.id) {				// Ensure that each processor runs once
			f.aopen() 						// Open for appending to file
			for i=0, pnm.idvec.size-1 {
				f.printf("%.3f %d\n", pnm.spikevec.x[i], pnm.idvec.x[i])	// Print the spike time and spiking cell gid
			}
			f.close()
		}
		pc.barrier()
	}
}

// also write out vector recordings
objref f, strobj
strobj = new StringFunctions()
strdef outfile, beforestr, afterstr
objref tgt

proc printtrace() {
	{sprint(outfile, "./netclamp/%s/%s/%s__%d_%s_trace_%03.0f.dat", BasedOnRunName, RunName, CellOI, gid, secstr, simnum)}
	{f = new File(outfile)}
	{f.wopen()}
	{f.printf("t\tv")} // write header
	if (myi_flag==1) {
		
		mt = new MechanismType(0)
		objref mechstring[mt.count()]
		k = 0
		for i=0, mt.count()-1 {
			mt.select( i )
			mt.selected(mname)
			if( ismembrane(mname)) {
				if (strcmp(mname,"capacitance")!=0 && strcmp(mname,"morphology")!=0 && strcmp(mname,"iconc_Ca")!=0  && strcmp(mname,"iconcCa")!=0 && strcmp(mname,"ccanl")!=0 && strcmp(mname,"cad")!=0 && strcmp(mname,"pas")!=0 ) { //
					if (strobj.substr(mname,"_ion")==-1) {
						//printf("myi_%s \n", mname) 
						sprint(tmpstr, "myi_%s", mname) // "cell.soma.%s(0.5)", tmpstr
						mechstring[k] = new String()
						mechstring[k].s = tmpstr
			if (strcmp("myi_ch_Nav", mechstring[k].s)==0) {
						print nickname, ": ", tmpstr}
						{k = k+1}
					} 
				} 
			}
		}
		mechlength.x[0]=k
		
		
		for rt = 0, mechlength.x[0]-1 {
			sprint(cmdstr, "%s_%s", nickname, mechstring[rt].s)		
		if (strcmp("myi_ch_Nav", mechstring[rt].s)==0) {
			print "Nav nickname: ", nickname}
		if (strcmp("som", nickname)==0) {
			print "cmdstr: ", cmdstr, " declared: ", name_declared(cmdstr)		}
			if (name_declared(cmdstr)>0) {
				f.printf("\t%s", mechstring[rt].s)
			}
		}
	}
	{f.printf("\n")}

	for i=0, (tstop-tstart)/dt-1 {
		{sprint(cmdstr, "{f.printf(\"%%g\\t%%g\", i*dt, mytrace%s.x[i])}", nickname)}
		{execute1(cmdstr)}
		if (myi_flag==1) {
			for rt = 0, mechlength.x[0]-1 {
				sprint(cmdstr, "%s_%s", nickname, mechstring[rt].s)				
				if (name_declared(cmdstr)>0) {
					sprint(cmdstr, "{f.printf(\"\\t%%g\", %s_%s.x[i])}", nickname, mechstring[rt].s)
					execute1(cmdstr)
				}
			}
		}						
		{f.printf("\n")}
	}
	{f.close()}	
}

use_cache_efficient=1
get_spike_hist=0
use_bin_queue=0
use_spike_compression=0
if (use_spike_compression==1) {
	maxstepval = 2.5
} else {
	maxstepval = 0.5
}														
															
proc rrun(){									// Run the network simulation and write out the results

	pnm.want_all_spikes() 						// Record all spikes of all cells on this machine into the
												//  vectors pnm.spikevec (spiketimes) and pnm.idvec (gids)
												
	local_minimum_delay = pc.set_maxstep(maxstepval)	// Set every machine's max step size to minimum delay of
												//  all netcons created on pc using pc.gid_connect, but
												//  not larger than 10
												
	stdinit()									// Call the init fcn (which is redefined in this code) and
												//  then make other standard calls (to stop watch, etc)

	pc.psolve(tstop)							// Equivalent to calling cvode.solve(tstop) but for parallel NEURON;
												//  solve will be broken into steps determined by the result of
												//  set_maxstep

	spikeout()	// Write the spike times (in ms) and spiking cells (by gid) into a file called "spikeraster.dat"
	
	printtrace()
}


walltime = startsw()
strdef cmd, cmdo
newtstop = tstop
warningflag=0

{cvode.cache_efficient(use_cache_efficient)} // always double check that this addition does not affect the spikeraster (via pointers in mod files, etc)

if (use_bin_queue==1) {
	use_fixed_step_bin_queue = 1 // boolean
	use_self_queue = 0 // boolean - this one may not be helpful for me, i think it's best for large numbers of artificial cells that receive large numbers of inputs
	{mode = cvode.queue_mode(use_fixed_step_bin_queue, use_self_queue)}
}

if (maxstepval==1) {
	nspike = 3 // compress spiketimes or not
	gid_compress = 0 //only works if fewer than 256 cells on each proc
	{nspike = pc.spike_compress(nspike, gid_compress)}
}
celsius = 34

simnum=0
numsims=10
lowerbound=0 // 5
rrun()

strdef f2n, f3n
objref f2, f3
f2 = new File()
f3 = new File()
sprint(f2n, "./netclamp/%s/%s_%d_NetStimConns.dat", BasedOnRunName, CellOI,  gid)
sprint(f3n, "./netclamp/%s/%s_%d_key.dat", BasedOnRunName, CellOI,  gid)
f3.wopen(f3n)
for simnum=1, numsims {
	//nclist.remove_all
	pregid=1

	f2.ropen(f2n)
	numStims = f2.scanvar

	numConnTypes = f2.scanvar
	postcellType = CellOItype
	
	f3.printf("%d\t%f\n", simnum, 1 + (simnum-lowerbound)*.1)

	for r=0, numConnTypes-1 {
		precellType = f2.scanvar
		numSynTypes = f2.scanvar
		synWeight = cellType[precellType].wgtConns.x[postcellType]	
		
		if (strcmp(cellType[precellType].cellType_string,"ca3cell")==0 || strcmp(cellType[precellType].cellType_string,"eccell")==0) {
			synWeight = synWeight*(1 + (simnum-lowerbound)*.1)
		}
			
		for SynNumber = 0, numSynTypes-1 {		
			nclist.o(pregid-1).weight = synWeight
			print "pregid = ", pregid, ", precellType = ", precellType, ", synapse = ", SynNumber, ", weight = ", synWeight
			//nc_append(pregid, cellType[postcellType].cellStartGid, precellType, SynNumber, synWeight + (pregid+1)*1000, 0.5)	// Make the connection
			pregid += 1
			cellType[postcellType].numCons.x[precellType] += 1
		}
	}

	f2.close()

	rrun()	// Run the network simulation (in proc rrun)
}

f3.close()

///////////////

{pc.runworker()} 	// Everything after this line is executed only by the host node
					//  The NEURON website describes this as "The master host returns immediately. Worker hosts
					//  start an infinite loop of requesting tasks for execution." 
					
{pc.done()}			// Sends the quit message to the worker processors, which then quit NEURON

quit()	// Sends the quit message to the host processor, which then quits NEURON

Loading data, please wait...