Getting the full topology (segments level) from hoc

Anything that doesn't fit elsewhere.
roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Getting the full topology (segments level) from hoc

Post by roybens » Fri Nov 08, 2013 1:02 am

Hi,
Is there a way to get the whole topology as in which segment compartments are connected?
I saw that "nt._v_parent_index" holds this information but how can i get to it from the hoc file

Another question:
I am trying to make a modify the dendritic tree so there are no segments that has more than one child. So i am counting the number of children from the original tree and changing the parents number of segments accordingly. Then i am disconnecting all children and reconnect them at 1/nsegs intervals. Though when i look at _v_parent_index there are still sections that are connected to the same segment in the original parent what am i doing wrong?.
my function looks like this:

Code: Select all

proc disperseBranching(){localobj sl,conVec,tempSec,existList strdef fatherStr
	forall{
		comps = 0
		sl=new SectionList()
		existList =new SectionList()
		sl.children()
		somaref= new SectionRef()
		fatherStr = secname()
		forsec sl{
			comps+=1
		}
		if (comps>somaref.sec.nseg){
			somaref.sec.nseg = comps
		}
		dseg = 1/somaref.sec.nseg
		connCount=somaref.sec.nseg
		forsec sl{
			newCon = dseg*connCount
			connCount-=1
			printf("newCon is %f  dseg is %f currently in %s\n",newCon,dseg,secname())
			disconnect()
			secref= new SectionRef()
			somaref.sec connect secref.sec(0), (newCon)
		}	
	printf("finished with %s\n",fatherStr)
	}
	
	
}
Thanks!

ted
Site Admin
Posts: 5351
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the full topology (segments level) from hoc

Post by ted » Fri Nov 08, 2013 11:45 pm

roybens wrote:Is there a way to get the whole topology as in which segment compartments are connected?
I saw that "nt._v_parent_index" holds this information but how can i get to it from the hoc file
I really don't understand what you are asking. First, I have no idea what nt._v_parent_index is supposed to be. Second, and actually more important: the recommended practice is that the 0 end of the child section should be connected to the 1 end of its parent section, except when the parent section is the root section (the section that has no parent), in which case the child may be connected to the parent's 0 end. A child section should never be connected to a parent's internal node, because the number and locations of internal nodes in a section depend on that section's nseg parameter, and one must be free to select any value for nseg in order to achieve good spatial resolution, without having to be concerned about what the value of nseg might do to where child branches are attached.

For this reason the CellBuilder cannot be used to create a model in which a child branch is attached to an internal node of its parent, and the Import3D tool always uses branch points to define the 0 and 1 ends of sections.

Of course, by writing one's own code it is always possible to do all kinds of things, even things that violate recommended practice. But if you do that, you're on your own.
I am trying to make a modify the dendritic tree so there are no segments that has more than one child.
Not necessary if you follow the recommended practice. And if you follow the recommended practice, then by definition there can be no branch point at which the parent "segment" (the segment at the 1 end of the parent section) does not have more than one child. "Oh, but what about the 0 end of the root section? Can't that have just one child?" Yes it can have just one child, but then it's not really a branch point, is it?

If this does not set things straight, please explain what your goal is and maybe I'll be able to give you a more useful answer.

roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Re: Getting the full topology (segments level) from hoc

Post by roybens » Sat Nov 09, 2013 2:00 pm

Hi Ted,
Thanks for your time,
I probably should have been more clear regarding nt._v_parent_index so its a member of the struct NrnThread defined in nrnoc/multicore.h that holds an array of all the indices of the parents for each segment.
My goal is to use GPU and CUDA to simulate neuronal models (http://www.frontiersin.org/Neuroinforma ... 4/abstract). So as our implementation in CUDA would be much faster if each segment has not more than two children. So i understand that it is not the recommended practice but still i am interested in finding a way to manipulate the tree so segments would have max two children. So using the procedure above should have brought this result but it did'nt.

Thanks
Roy

hines
Site Admin
Posts: 1523
Joined: Wed May 18, 2005 3:32 pm

Re: Getting the full topology (segments level) from hoc

Post by hines » Thu Nov 14, 2013 10:10 am

The only way I can imagine doing that is to insert zero area and resistance nodes whereever there are more than two connections and then reconnect two child nodes to the new node. For example
if the soma has 5 dendrites, you would need a chain of 4 additional nodes. As long as you do this following the constraint that the parent index < node index, everyting should work properly.

roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Re: Getting the full topology (segments level) from hoc

Post by roybens » Fri Nov 15, 2013 12:44 pm

Thank you very much for this,

I tried to implement this but when comparing it to the original code (without splitting the branching) i get different results. So i probably got something wrong:
So when i set Ra=0 in the additional nodes then the volt is nan, so i changed it to 0.00000001 but the results between the voltages are very different.
I am pasting three hoc files split.hoc where the split procedure is applied, no_split.hoc where the code for splitting is commented out ,and sp.hoc which is a simple morphology file.
Thanks a lot

no_split.hoc

Code: Select all

objref somaref,secref,extVec,strFunc,stim,g1
strFunc = new StringFunctions()
load_file("./sp.hoc")
access soma
stim = new IClamp(0.5)
stim.del=10
stim.dur=100
stim.amp = 1
finitialize(-65)
tstop=300
g1 = new Graph()
g1.beginline()
forsec "soma"{
insert hh
gnabar_hh = 2*gnabar_hh
}
forsec "apic"{
insert hh
}
forsec "dend"{
insert hh
}
while(t<tstop){
	fadvance()
	g1.line(t,v(0.5))
}

split.hoc

Code: Select all

objref somaref,secref,extVec,strFunc,stim,g1
strFunc = new StringFunctions()
create Exts[1]
proc disperseBranch(){localobj sl,conVec,tempSec,existList strdef fatherStr
	ext_ind=0
	forall{
		comps = 0
		fatherStr = secname()
		if (strFunc.len(fatherStr)<4 || strFunc.substr(fatherStr,"Exts")==-1){//we dont want to iterate over Exts compartments 
			sl=new SectionList()
			sl.children()
			somaref= new SectionRef()
			print "priniting Children"
			forsec sl{
				comps+=1
			}
			if (comps>1){
				//First connecting between the extension (if we have 3 comps we need two connections)
				for(i=ext_ind;i<ext_ind+comps-1;i+=1){
				 connect Exts[i+1](0), Exts[i](1)
					 Exts[i]{
					L=0.00000001
					diam = 0.000001
					Ra = 0
					}
				}
				printf("printing ext_ind %d\n",ext_ind)
				somaref.sec connect Exts[ext_ind](0), 1
				forsec sl{
					disconnect()
					tempSec = new SectionRef()
					connect tempSec.sec(0), Exts[ext_ind](1)
					ext_ind+=1
				}
			}
		}
	}
}
proc countExts(){localobj sl
extVec = new Vector()
forall{
	comps = 0
	sl=new SectionList()
	sl.children()
	somaref= new SectionRef()
	forsec sl{
		comps+=1
	}
	if (comps>1){
		extVec.append(comps)
	}
	}
}
proc printChildren(){localobj sl
	sl = new SectionList()
	sl.children()
	u=0
	forsec sl{
	print secname()
	u+=1
	}
	print u
}	
proc printparent(){localobj sref
	sref = new SectionRef()
	sref.parent
	print secname()
}		
proc countSegs(){
segs=0
comps=0
	forall{
		segs+=nseg
		comps+=1
	}
}
load_file("./sp.hoc")
countExts()
ext_num = extVec.sum()
printf("Creating Exts %d\n",ext_num)
create Exts[ext_num]
disperseBranch()
access soma
stim = new IClamp(0.5)
stim.del=10
stim.dur=100
stim.amp = 1
finitialize(-65)
tstop=300
g1 = new Graph()
g1.beginline()
forsec "soma"{
insert hh
gnabar_hh = 2*gnabar_hh
}
forsec "apic"{
insert hh
}
forsec "dend"{
insert hh
}
while(t<tstop){
	fadvance()
	g1.line(t,v(0.5))
}
sp.hoc

Code: Select all

proc celldef() {
  topol()
  subsets()
  geom()
  biophys()
  geom_nseg()
}

create soma, dend[2], apic[2]

proc topol() { local i
  for i = 0, 1 connect dend[i](0), soma(0.5)
  connect apic(0), soma(0.5)
 connect apic[1](0), apic(1)
  basic_shape()
}
proc shape3d_1() {
  soma {pt3dclear()
	pt3dadd(-11.0933, -8.95771, 0, 1.83129)
	pt3dadd(-10.7893, -8.16715, 0, 3.07228)
	pt3dadd(-10.4853, -7.37658, 0, 3.7064)
	pt3dadd(-10.1814, -6.58602, 0, 5.16232)
	pt3dadd(-9.87739, -5.79545, 0, 7.5116)
	pt3dadd(-9.57342, -5.00488, 0, 8.60502)
	pt3dadd(-9.26945, -4.21432, 0, 9.645)
	pt3dadd(-8.96548, -3.42375, 0, 11.1731)
	pt3dadd(-8.66151, -2.63319, 0, 12.2908)
	pt3dadd(-8.35754, -1.84262, 0, 13.3186)
	pt3dadd(-8.05357, -1.05205, 0, 14.3568)
	pt3dadd(-7.7496, -0.261487, 0, 15.4776)
	pt3dadd(-7.44563, 0.529079, 0, 16.3335)
	pt3dadd(-7.14166, 1.31964, 0, 16.0558)
	pt3dadd(-6.83769, 2.11021, 0, 15.2669)
	pt3dadd(-6.53372, 2.90078, 0, 14.306)
	pt3dadd(-6.22975, 3.69134, 0, 11.799)
	pt3dadd(-5.92578, 4.48191, 0, 10.2575)
	pt3dadd(-5.62181, 5.27248, 0, 8.97583)
	pt3dadd(-5.31784, 6.06304, 0, 6.01391)
	pt3dadd(-5.01387, 6.85361, 0, 3.51222)
  }
  dend {pt3dclear()
	pt3dstyle(1, -7.82441, -0.45605, 0)
	pt3dadd(-5.61, 6.42, 0, 0.4)
	pt3dadd(-6.42, 8.83, 0, 0.4)
	pt3dadd(-7.62, 10.43, 0, 0.4)
	pt3dadd(-8.02, 12.44, 0, 0.4)
	pt3dadd(-8.42, 15.24, 0, 0.4)
	pt3dadd(-8.82, 16.85, 0, 0.4)
	pt3dadd(-8.82, 17.25, 0, 0.4)
	pt3dadd(-9.63, 19.25, 0, 0.4)
	pt3dadd(-9.63, 20.46, 0, 0.4)
	pt3dadd(-9.63, 20.86, 0, 0.4)
	pt3dadd(-10.03, 23.27, 0, 0.4)
	pt3dadd(-10.43, 26.48, 0, 0.4)
  }
  dend[1] {pt3dclear()
	pt3dstyle(1, -7.82441, -0.45605, 0)
	pt3dadd(-9.63, 8.42, 0, 0.4)
	pt3dadd(-9.63, 8.83, 0, 0.4)
	pt3dadd(-12.03, 12.44, 0, 0.4)
	pt3dadd(-12.83, 12.44, 0, 0.4)
	pt3dadd(-13.64, 12.84, 0, 0.4)
	pt3dadd(-16.84, 15.64, 0, 0.4)
	pt3dadd(-17.25, 16.05, 0, 0.4)
	pt3dadd(-18.45, 17.25, 0, 0.4)
	pt3dadd(-18.85, 17.65, 0, 0.4)
	pt3dadd(-20.86, 20.06, 0, 0.4)
	pt3dadd(-21.66, 20.46, 0, 0.4)
	pt3dadd(-22.46, 22.46, 0, 0.4)
	pt3dadd(-22.86, 23.27, 0, 0.4)
	pt3dadd(-23.26, 24.87, 0, 0.4)
	pt3dadd(-24.06, 26.88, 0, 0.4)
	pt3dadd(-24.06, 27.28, 0, 0.4)
  }
  apic {pt3dclear()
	pt3dstyle(1, -7.82441, -0.45605, 0)
	pt3dadd(1.2, -1.2, 0.15, 1.6)
	pt3dadd(1.6, -1.2, 0.15, 1.6)
	pt3dadd(2.81, -1.2, 0.15, 1.6)
	pt3dadd(3.61, -1.2, 0.15, 1.6)
	pt3dadd(4.81, -1.2, 0.15, 1.6)
	pt3dadd(5.21, -1.2, 0.15, 1.6)
	pt3dadd(6.42, -1.2, 0.15, 1.6)
	pt3dadd(8.82, 0, 0.15, 1.6)
	pt3dadd(9.22, 0, 0.15, 1.6)
	pt3dadd(11.23, 1.6, 0.15, 1.2)
	pt3dadd(11.23, 3.61, 0.15, 1.2)
	pt3dadd(11.63, 3.61, 0.15, 1.2)
	pt3dadd(13.24, 1.6, 0.15, 1.2)
	pt3dadd(13.64, 1.2, 0.15, 1.2)
	pt3dadd(16.84, 0, 0.15, 1.2)
	pt3dadd(17.25, -0.8, 0.15, 1.2)
	pt3dadd(17.65, -0.8, 0.15, 1.2)
	pt3dadd(20.45, -0.8, 0.15, 1.2)
	pt3dadd(20.86, -0.8, 0.15, 1.2)
	pt3dadd(24.47, -0.4, 0.15, 1.2)
	pt3dadd(25.27, -0.4, 0.15, 1.2)
	pt3dadd(27.27, -0.4, 0.15, 1.2)
	pt3dadd(27.67, -0.4, 0.15, 1.2)
	pt3dadd(28.88, 0.8, -12.21, 0.8)
	pt3dadd(30.48, -0.8, -12.21, 0.8)
	pt3dadd(33.29, -3.21, -12.21, 0.8)
	pt3dadd(33.69, -3.61, -12.21, 0.8)
	pt3dadd(35.29, -4.01, -12.21, 0.8)
	pt3dadd(35.7, -4.41, -12.21, 0.8)
	pt3dadd(37.3, -4.41, -12.21, 0.8)
	pt3dadd(38.1, -4.41, -12.21, 0.8)
	pt3dadd(40.91, -4.01, -12.21, 0.8)
	pt3dadd(41.31, -4.01, -12.21, 0.8)
	pt3dadd(43.72, -2.81, -12.21, 0.8)
	pt3dadd(44.12, -2.81, -12.21, 0.8)
	pt3dadd(46.52, -2.01, -12.21, 0.8)
	pt3dadd(47.33, -1.6, -12.21, 0.4)
	pt3dadd(47.73, -1.6, -12.21, 0.4)
	pt3dadd(50.13, -1.6, -12.21, 0.4)
	pt3dadd(50.53, -1.6, -12.21, 0.4)
	pt3dadd(53.34, -0.4, -12.21, 0.4)
	pt3dadd(54.14, -0.4, -12.21, 0.4)
	pt3dadd(56.15, 0.4, -12.21, 0.4)
	pt3dadd(56.55, 0.4, -12.21, 0.4)
	pt3dadd(58.16, 0.4, -12.21, 0.4)
	pt3dadd(60.96, -1.6, -12.21, 0.4)
	pt3dadd(62.17, -2.41, -12.21, 0.4)
	pt3dadd(62.57, -2.41, -12.21, 0.4)
	pt3dadd(64.57, -2.41, -12.21, 0.4)
	pt3dadd(65.78, -2.41, -12.21, 0.4)
	pt3dadd(66.18, -2.01, -12.21, 0.4)
	pt3dadd(67.38, -1.2, -12.21, 0.4)
	pt3dadd(67.78, -1.2, -12.21, 0.4)
	pt3dadd(68.98, -0.4, -12.21, 0.4)
	pt3dadd(69.79, -0.4, -12.21, 0.4)
	pt3dadd(72.99, 0, -12.21, 0.4)
	pt3dadd(73.8, 0, -12.21, 0.4)
	pt3dadd(76.2, 0, -12.21, 0.4)
	pt3dadd(76.6, 0, -12.21, 0.4)
	pt3dadd(78.61, 0, -12.21, 0.4)
	pt3dadd(79.41, 0, -12.21, 0.4)
	pt3dadd(83.02, 0, -12.21, 0.4)
	pt3dadd(83.42, 0, -12.21, 0.4)
	pt3dadd(86.23, 0, -12.21, 0.4)
	pt3dadd(87.43, 1.2, -12.21, 0.8)
	pt3dadd(91.85, -0.4, -12.21, 0.4)
	pt3dadd(95.86, -0.4, -12.21, 0.4)
	pt3dadd(96.26, -0.4, -12.21, 0.4)
	pt3dadd(99.06, 1.2, -12.21, 0.4)
	pt3dadd(102.67, 2.01, -12.21, 0.4)
	pt3dadd(103.88, 2.81, -13.85, 5.61)
	pt3dadd(106.28, 4.41, -13.85, 1.6)
	pt3dadd(108.29, 5.21, -17.09, 1.2)
	pt3dadd(110.29, 5.21, -17.09, 1.2)
	pt3dadd(111.1, 5.21, -17.09, 1.2)
	pt3dadd(113.9, 4.81, -17.09, 1.2)
	pt3dadd(114.71, 4.41, -17.09, 1.2)
  }
  apic[1] {pt3dclear()
	pt3dadd(114.71, 4.41, -17.09, 0.4)
	pt3dadd(130.75, 3.61, -20.73, 0.4)
	pt3dadd(141.18, 2.81, -20.73, 0.4)
	pt3dadd(147.99, 3.61, -20.73, 0.4)
	pt3dadd(162.83, 3.21, -20.73, 0.4)
	pt3dadd(170.45, 2.81, -20.73, 0.4)
	pt3dadd(170.86, 2.81, -20.73, 0.4)
	pt3dadd(178.48, 3.61, -20.73, 0.4)
	pt3dadd(186.1, 4.01, -20.73, 0.4)
	pt3dadd(186.5, 4.01, -20.73, 0.4)
	pt3dadd(190.91, 6.02, -20.73, 0.4)
	pt3dadd(191.31, 6.02, -20.73, 0.4)
	pt3dadd(195.72, 7.22, -20.73, 0.4)
	pt3dadd(196.12, 7.22, -20.73, 0.4)
	pt3dadd(201.34, 6.82, -20.73, 0.4)
	pt3dadd(201.74, 6.82, -20.73, 0.4)
	pt3dadd(206.95, 6.02, -20.73, 0.4)
	pt3dadd(207.35, 6.02, -20.73, 0.4)
	pt3dadd(211.76, 5.62, -20.73, 0.4)
	pt3dadd(221.79, 4.81, -20.73, 0.4)
  }
 
}
proc basic_shape() {
  shape3d_1()
}

objref all, somatic, basal, apical
proc subsets() { local i
  objref all, somatic, basal, apical
  all = new SectionList()
    soma all.append()
    for i=0, 1 dend[i] all.append()
    for i=0, 1 apic[i] all.append()

  somatic = new SectionList()
    soma somatic.append()

  basal = new SectionList()
    for i=0, 1 dend[i] basal.append()

  apical = new SectionList()
    for i=0, 1 apic[i] apical.append()

}
proc geom() {
}
proc geom_nseg() {
}
proc biophys() {
}
access soma

celldef()

ted
Site Admin
Posts: 5351
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the full topology (segments level) from hoc

Post by ted » Sun Nov 17, 2013 12:08 pm

Does this help?

Code: Select all

oc>topology()

|-|       soma(0-1)
   `|       Exts[0](0-1)
     `|       Exts[1](0-1)
       `|       Exts[2](0-1)
         `|       dend[0](0-1)
       `|       dend[1](0-1)
     `|       apic[0](0-1)
       `|       apic[1](0-1)

oc>for ii=0,2 print Exts[ii].L
first instance of ii
1e-08 
1e-08 
100 
oc>forsec "Exts" print diam
1e-06 
1e-06 
500 

roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Re: Getting the full topology (segments level) from hoc

Post by roybens » Sun Nov 17, 2013 2:07 pm

Yes! it does help a lot, my bad
Thanks so much!

ted
Site Admin
Posts: 5351
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the full topology (segments level) from hoc

Post by ted » Sun Nov 17, 2013 2:27 pm

Good that you were using a "toy model," which makes testing easier.

When you have worked it out, what's the fix?

roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Re: Getting the full topology (segments level) from hoc

Post by roybens » Sun Nov 17, 2013 4:49 pm

So in the loop i didn't change the properties of the last Ext section i am pasting the corrected code:

Code: Select all




objref somaref,secref,extVec,strFunc,stim,g1
strFunc = new StringFunctions()
create Exts[1]
proc disperseBranch(){localobj sl,conVec,tempSec,existList strdef fatherStr
	ext_ind=0
	forall{
		comps = 0
		fatherStr = secname()
		if (strFunc.len(fatherStr)<4 || strFunc.substr(fatherStr,"Exts")==-1){//we dont want to iterate over Exts compartments 
			sl=new SectionList()
			sl.children()
			somaref= new SectionRef()
			print "priniting Children"
			forsec sl{
				comps+=1
			}
			if (comps>1){
				//First connecting between the extension (if we have 3 comps we need two connections)
				for(i=ext_ind;i<ext_ind+comps;i+=1){
				if(i<ext_ind+comps-1){
				 connect Exts[i+1](0), Exts[i](1)
				 }
					 Exts[i]{
					L=0.000000000000001
					diam = 0.0000000000001
					Ra = 0.000000000000001
					}
				}
				printf("printing ext_ind %d\n",ext_ind)
				somaref.sec connect Exts[ext_ind](0), 1
				forsec sl{
					disconnect()
					tempSec = new SectionRef()
					connect tempSec.sec(0), Exts[ext_ind](1)
					ext_ind+=1
				}
			}
		}
	}
}
proc countExts(){localobj sl
extVec = new Vector()
forall{
	comps = 0
	sl=new SectionList()
	sl.children()
	somaref= new SectionRef()
	forsec sl{
		comps+=1
	}
	if (comps>1){
		extVec.append(comps)
	}
	}
}
proc printChildren(){localobj sl
	sl = new SectionList()
	sl.children()
	u=0
	forsec sl{
	print secname()
	u+=1
	}
	print u
}	
proc printparent(){localobj sref
	sref = new SectionRef()
	sref.parent
	print secname()
}		
proc countSegs(){
segs=0
comps=0
	forall{
		segs+=nseg
		comps+=1
	}
}
load_file("./sp.hoc")
countExts()
ext_num = extVec.sum()
printf("Creating Exts %d\n",ext_num)
create Exts[ext_num]
disperseBranch()
access soma
stim = new IClamp(0.5)
stim.del=10
stim.dur=100
stim.amp = 1
finitialize(-65)
tstop=300
g1 = new Graph()
g1.beginline()
forsec "soma"{
insert hh
gnabar_hh = 2*gnabar_hh
}
forsec "apic"{
insert hh
}
forsec "dend"{
insert hh
}
while(t<tstop){
	fadvance()
	g1.line(t,v(0.5))
}
This solution is great, Though there is still a ~0.1 mv error between the two simulation when i am using the toy model, in actual model (Mainen & Sejnowski (1996)) i get a larger error of ~3.5 mv
It arises probably since you cannot actually have 0 area and resistance section. I could live with this error, though if you guys have an idea how i can make it better it would be great.
Thanks a lot for this solution
Roy

ted
Site Admin
Posts: 5351
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Getting the full topology (segments level) from hoc

Post by ted » Sun Nov 17, 2013 5:18 pm

Not a good idea to make the diameters small--multiplying diam by 1/k is equivalent to multiplying Ra by k^2. Also, don't forget about cm. Suggest
L = 1e-9
diam = 1
Ra = 1e-9
cm = 1e-9

hines
Site Admin
Posts: 1523
Joined: Wed May 18, 2005 3:32 pm

Re: Getting the full topology (segments level) from hoc

Post by hines » Sun Nov 17, 2013 5:37 pm

You should be able to get quantitative identiy after your transformation if you set the connection coefficients for your extra nodes using an approach analogous to that of
https://senselab.med.yale.edu/modeldb/S ... mod\ri.mod
the only thing of importance is to set
NODEA(nd) = ...;
NODEB(nd) *= ...;
properly.
There are two considerations. 1) All the introduced nodes have equations which end up so that the resullt is vchild = vparent
2) all the original child NODEA and NODEB will have the same values if they were originally connected to a 0 or 1 end of the parent.

One minor administrative suggestion. If you have n children, You only need to disconnect and insert a single ext section with nseg = n -1 (A section always has nseg + 1 node where the last is
a 0 area node) Then connect one child to the nseg+1 locations of the ext section. eg.
i=0
ext for (x) (if x > 0) { connect child(0), x i += 1 }
Then call the function in your connection_coef.mod file which sets the NODEA and NODEB of each segment in the ext section. There may be a few other things to do on the original parent and
child nodes but I don't think so as long as you consider that those are current balance equations each divided by the area of the node (except for zero area nodes). Ie you must formulate
each of your ext node equations as a sum of three currents = 0

hines
Site Admin
Posts: 1523
Joined: Wed May 18, 2005 3:32 pm

Re: Getting the full topology (segments level) from hoc

Post by hines » Sun Nov 17, 2013 5:38 pm

sorry. meant 'NODEB = ...'

roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Re: Getting the full topology (segments level) from hoc

Post by roybens » Mon Nov 18, 2013 5:59 pm

regarding the verbatim used in the ri.mod it helps me a lot since now i can get the fmatrix() in any precision i like and the segments level topology (_nt._v_parent_index) so it has been super helpful and answers the question i asked in the topic.
here is the mod file i am using now:

Code: Select all


NEURON { SUFFIX nothing }

VERBATIM
char* secname();
ENDVERBATIM

PROCEDURE scale_connection_coef(x, factor) {
VERBATIM {
	Section* sec;
	Node* nd;
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif

	sec = chk_access();
	if (_lx <= 0. || _lx > 1.) {
		hoc_execerror("out of range, must be 0 < x <= 1", (char*)0);
	}
	/*printf("scale_connection_coefs %s(%g) %d\n", secname(sec), _lx, sec->nnode);*/
	/* assumes 0 end of child connected to parent */
	if (_lx == 1.) {
		nd = sec->pnode[sec->nnode-1];
	}else{
		nd = sec->pnode[(int) (_lx*(double)(sec->nnode-1))];
	}
	/*printf("%g %g\n", NODEA(nd), NODEB(nd));*/
#if defined(t)
	_nt = nd->_nt;
#endif
	NODEA(nd) *= _lfactor;
	NODEB(nd) *= _lfactor;
}
ENDVERBATIM
}
PROCEDURE MyPrintMatrix() {
VERBATIM {
	Section* sec;
	Node* nd;
	int ii;
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif
for(ii=0;ii<_nt->end;ii++){
nd=_nt->_v_node[ii];
printf("%d %1.15f %1.15f %1.15f %1.15f\n", ii, NODEB(nd), NODEA(nd), NODED(nd), NODERHS(nd));
}

}
ENDVERBATIM
}
PROCEDURE MyTopology() {
VERBATIM {
	int ii;
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif
for(ii=0;ii<_nt->end;ii++){

printf("%d %d\n", ii, _nt->_v_parent_index[ii]);
}
}
ENDVERBATIM
}
I am still having some problem regarding the issue of reconnecting the tree so each segment would have no more than two children:
I want to implement the solution you suggested connecting the n children to a section with n-1 segments and then do the corrections in the relevant NODEA and NODEB.
But when i am connecting somas' children i get an error:
nrniv: dend[0] must connect at position 0 or 1
so i probably did something wrong.
Any ideas? (I am adding the updated code for split.hoc)
Anyhow Thanks so much for your time this is being extremely helpful

split.hoc

Code: Select all

objref somaref,secref,extVec,strFunc,stim,g1
strFunc = new StringFunctions()
create Exts[1]
proc hinesDisperseBranching(){localobj sl,secRef
	counter = 0
	forall{
		secRef = new SectionRef()
		if(secRef.nchild>1){
			Exts[counter]{
				L = 1e-9
				diam = 1
				Ra = 1e-9
				cm = 1e-9
				nseg = secRef.nchild-1
			}
			i=0
			Exts[counter]{
				for(x){ 
					if (x > 0) {
						printf("x is %f i is %d\n",x,i)
						connect secRef.child[i](0), x
						i +=1
					}
				}
			}
			sl = new SectionList()
			sl.children()
			forsec sl{
				disconnect()
			}
	counter+=1		
		}
	}
}
proc countExts(){
ext_num=0
forall{
	somaref= new SectionRef()
	if (somaref.nchild>1){
		ext_num+=1
	}
	}
}
load_file("./sp.hoc")
countExts()
printf("Creating Exts %d\n",ext_num)
create Exts[ext_num]
hinesDisperseBranching()
access soma
stim = new IClamp(0.5)
stim.del=10
stim.dur=100
stim.amp = 1
finitialize(-65)
tstop=300
g1 = new Graph()
g1.beginline()
forsec "soma"{
insert hh
gnabar_hh = 2*gnabar_hh
}
forsec "apic"{
insert hh
}
forsec "dend"{
insert hh
}
while(t<tstop){
	fadvance()
	g1.line(t,v(0.5))
}
Roy

hines
Site Admin
Posts: 1523
Joined: Wed May 18, 2005 3:32 pm

Re: Getting the full topology (segments level) from hoc

Post by hines » Mon Nov 18, 2013 7:54 pm

shouldn't you check that the number of children is > 2 instead of > 1

roybens
Posts: 39
Joined: Fri Mar 14, 2008 7:57 am

Re: Getting the full topology (segments level) from hoc

Post by roybens » Mon Nov 18, 2013 10:26 pm

I am checking if nchild>1 so if the section does have more than one children there is a scenario where there are two children in the same segment(now that i can use the verbatim i should probably do a better check - that there is no segment that appears more than twice in the v_parent_node array).
But even if i change the condition to
if(secRef.nchild>2)
i still get the same error at the soma after connecting three out of the four children to Exts[0]

Code: Select all

x is 1.000000 i is 2
nrniv: dend[0] must connect at position 0 or 1

Post Reply