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 }