A single column thalamocortical network model (Traub et al 2005)

 Download zip file   Auto-launch 
Help downloading and running models
To better understand population phenomena in thalamocortical neuronal ensembles, we have constructed a preliminary network model with 3,560 multicompartment neurons (containing soma, branching dendrites, and a portion of axon). Types of neurons included superficial pyramids (with regular spiking [RS] and fast rhythmic bursting [FRB] firing behaviors); RS spiny stellates; fast spiking (FS) interneurons, with basket-type and axoaxonic types of connectivity, and located in superficial and deep cortical layers; low threshold spiking (LTS) interneurons, that contacted principal cell dendrites; deep pyramids, that could have RS or intrinsic bursting (IB) firing behaviors, and endowed either with non-tufted apical dendrites or with long tufted apical dendrites; thalamocortical relay (TCR) cells; and nucleus reticularis (nRT) cells. To the extent possible, both electrophysiology and synaptic connectivity were based on published data, although many arbitrary choices were necessary.
1 . Traub RD, Contreras D, Cunningham MO, Murray H, LeBeau FE, Roopun A, Bibbig A, Wilent WB, Higley MJ, Whittington MA (2005) Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. J Neurophysiol 93:2194-232 [PubMed]
2 . Traub RD, Contreras D, Whittington MA (2005) Combined experimental/simulation studies of cellular and network mechanisms of epileptogenesis in vitro and in vivo. J Clin Neurophysiol 22:330-42 [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: Neocortex; Thalamus;
Cell Type(s): Thalamus geniculate nucleus/lateral principal GLU cell; Thalamus reticular nucleus GABA cell; Neocortex U1 L6 pyramidal corticalthalamic GLU cell; Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex fast spiking (FS) interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s): I Na,p; I Na,t; I L high threshold; I T low threshold; I A; I K; I M; I h; I K,Ca; I Calcium; I A, slow;
Gap Junctions: Gap junctions;
Receptor(s): GabaA; AMPA; NMDA;
Simulation Environment: NEURON; FORTRAN;
Model Concept(s): Activity Patterns; Bursting; Temporal Pattern Generation; Oscillations; Simplified Models; Epilepsy; Sleep; Spindles;
Implementer(s): Traub, Roger D [rtraub at us.ibm.com];
Search NeuronDB for information about:  Thalamus geniculate nucleus/lateral principal GLU cell; Thalamus reticular nucleus GABA cell; Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex U1 L6 pyramidal corticalthalamic GLU cell; GabaA; AMPA; NMDA; I Na,p; I Na,t; I L high threshold; I T low threshold; I A; I K; I M; I h; I K,Ca; I Calcium; I A, slow;
Files displayed below are from the implementation
balcomp.hoc *
defvar.hoc *
lbcreate.hoc *
mscreate.hoc *
parlib2.hoc *
traubcon.hoc *
traubcon_net.hoc *
if (name_declared("use_traubexact") == 0) {
	use_traubexact = 0
begintemplate TraubExactInfo
public tci, traub_parent
objref tci[14], traub_parent[14]
endtemplate TraubExactInfo
objref traubExactInfo
if (use_traubexact) {	// for network context traub exact connections
	traubExactInfo = new TraubExactInfo()

objref tci, traub_parent
strdef tstr
proc traub_connect() {local x, ip, ic  localobj child
	if (numarg() == 1) {
		tci = new Matrix($1, $1, 2)
		traub_parent = new Vector($1)
	if (use_traubexact) {
		if (traubExactInfo.tci[$o1.type] == nil) {
			traubExactInfo.tci[$o1.type] = tci
			traubExactInfo.traub_parent[$o1.type] = traub_parent
	if (numarg() == 5) {x = $5} else { x = 1 }
	ip = $2  ic = $3
	// which is the child and which is the parent?
	// assume a child is connected to the parent before it
	// makes a delta connection (to child connected to same parent)
	// if $3 is already connected assume the child is $2
	$o1.comp[$3] child = new SectionRef()
	if (child.has_parent) {
		$o1.comp[$2] child = new SectionRef()
		ip = $3  ic = $2
	// standard NEURON connection only if child not already connected
	if (!child.has_parent) {
		$o1.comp[ip] connect child.sec(0), x
		traub_parent.x[ic] = ip
	// save info for optional traub style connections
	tci.x[$2][$3] = $4
	tci.x[$3][$2] = $4
	// note that the above matrix has the topology of a tree unless
	// child to child connections form delta loops.

// exact() transforms the delta and pyramidal axial resistance branch points to Y's and X's.
// single_cell points at current cell object, .e.g TCR[0]
proc exact(){ 
	print "exact() has been called."
	execute("traubexact(single_cell, tci)")

proc traubexact() {
	//$o1 is the cell, $o2 is the connection matrix
	// assume the names of sections range from
	// $o1.comp[1] to $o1.comp[$o2.nrow-1]
	// the connection matrix

	// the traub connection rule is that a childs 0 end is connected
	// to the .5 location of the parent and the resistance is defined
	// by the connection matrix (uS) with one exception. Sometimes
	// pairs or triples of children are also connected to each other
	// and these children are moved to the 1 end of the parent and
	// the delta network is replaced by the equivalent wye network.

	traubexact_topology($o1, $o2, traub_parent)
	traubexact_coef($o1, $o2, traub_parent)

proc traubexact_topology() { local j, jc, ip, ic, gc, x
	// move all children to .5 location of parent except if
	// a group of children have a delta topology in which case move
	// them to the 1 location of the parent.
	for ic = 2, $o2.nrow-1 $o1.comp[ic] {
		// parent must be the least index
		ip = $o3.x[ic]
		x = .5
		for j = 0, $o2.sprowlen(ic)-1 {
			gc = $o2.spgetrowval(ic, j, &jc)
			// the question is whether jc is also connected
			// to ic's parent.
			if ($o2.getval(ip, jc)) { // yes
				x = 1
		if (x != parent_connection()) {
			$o1.comp[ip] connect $o1.comp[ic](0), x

proc traubexact_coef() { local i, j, id, ip, ic, gp, gc, jc, x localobj dtopol, dmat, dindx, wye
	// now that the topology is stable we can safely change the
	// connection coefficients so they reflect the desired axial
	// resistance
	dtopol = new Vector($o2.nrow)
	for ic = 2, $o2.nrow-1 $o1.comp[ic] {
		ip = $o3.x[ic]
		gp = $o2.getval(ip,ic)
		if (parent_connection() == 0.5) {
			rold = ri(0.5)
			scale_connection_coef(0.5, rold*gp)
		}else{ // this is delta-topology child.
			// so save info for further analysis
			dtopol.x[ic] = ip
	// examine the dtopol info. For each delta topology, construct a
	// conductance matrix and a section index vector (the parent is always
	// the first index) and then find the equivalent wye conductances and
	// use those to specify the connection coefficients
	for (id = dtopol.indwhere(">",0); id != -1; id = dtopol.indwhere(">",0)){
		// a delta exists and the parent is
		ip = dtopol.x[id]
		// note: error if ip != $o3.x[id]
		// so all the children are
		dindx = dtopol.c.indvwhere("==", ip)
		// prepare for next iteration over id
		for i=0, dindx.size-1 { dtopol.x[dindx.x[i]] = 0 }
		// and the complete delta index vector is
		dindx.insrt(0, ip)		
		// the delta conductance matrix
		dmat = new Matrix(dindx.size, dindx.size)
		if (dindx.size < 3) {
			execerror("not a delta", 0)
		for i = 0, dindx.size-1 {
			ic = dindx.x[i]
			for j = i, dindx.size-1 {
				jc = dindx.x[j]
				dmat.x[i][j] = $o2.getval(ic,jc)
				dmat.x[j][i] = dmat.x[i][j]
		// transform into an equivalent wye
		wye = delta2wye(dmat)
		// modify the connection coefficients
		for i = 0, dindx.size-1 $o1.comp[dindx.x[i]] {
			if (i == 0) { x = 1 } else { x = .5 }
			rold = ri(x)
			scale_connection_coef(x, rold/wye.x[i])

obfunc delta2wye() {localobj g, x, row0
	// $o1 is the conductance matrix. return vector is the equivalent
	// wye vector of resistances.
	// construct a full Gij*(Vi - Vj) = Ii matrix equation
	x = new Vector($o1.nrow)
	$o1.mulv(x, x)
	g = new Matrix($o1.nrow, $o1.ncol)
	g.setdiag(0, x)
	//presently indeterminate. save first row and replace by ground equation
	row0 = g.getrow(0)
	g.setrow(0, 0)
	g.x[0][0] = 1
	$o1.getdiag(0, x) // 1 to n-1 are the r0 + ri resistances and x.x[0]=1
	// need one more equation (not involving rp) and assume
	// numerical accuracy allows an arbitrary pair of nodes for
	// current injection and ground.
	// so put the row back, select node 1 as the ground point,
	// and inject into node 2.
	g.setrow(0, row0)
	g.setrow(1, 0)
	g.x[1][1] = 1
	//wasteful but...
	x.x[0] = $o1.x[2][2]
	// now the 0 element is r1+r2
	// set up the matrix equation to solve for ri
	g.setcol(0, 1)
	g.setdiag(0, 1)
	// but the first equation is r1 + r2 = x0
	g.x[0][0] = 0
	g.x[0][1] = 1
	g.x[0][2] = 1
	x = g.solv(x)

	return x