Cell splitting in neural networks extends strong scaling (Hines et al. 2008)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:97917
Neuron tree topology equations can be split into two subtrees and solved on different processors with no change in accuracy, stability, or computational effort; communication costs involve only sending and receiving two double precision values by each subtree at each time step. Application of the cell splitting method to two published network models exhibits good runtime scaling on twice as many processors as could be effectively used with whole-cell balancing.
Reference:
1 . Hines ML, Eichner H, Schürmann F (2008) Neuron splitting in compute-bound parallel network simulations enables runtime scaling with twice as many processors. J Comput Neurosci 25:203-10 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Generic;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Methods;
Implementer(s): Hines, Michael [Michael.Hines at Yale.edu];
/
splitcell
nrntraub
hoc
balcomp.hoc *
binfo.hoc *
defvar.hoc *
karkar.hoc
lbcreate.hoc *
loadbal.hoc *
mscreate.hoc
msdiv.hoc *
parlib.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)
		return
	}
	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)
	doNotify()
	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
				break
			}
		}
		if (x != parent_connection()) {
			disconnect()
			$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)
	x.fill(1)
	$o1.mulv(x, x)
	g = new Matrix($o1.nrow, $o1.ncol)
	g.setdiag(0, x)
	g.add($o1.muls(-1))
	//presently indeterminate. save first row and replace by ground equation
	row0 = g.getrow(0)
	g.setrow(0, 0)
	g.x[0][0] = 1
	g.inverse($o1)
	$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...
	g.inverse($o1)
	x.x[0] = $o1.x[2][2]
	// now the 0 element is r1+r2
	// set up the matrix equation to solve for ri
	g.zero
	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
}

Loading data, please wait...