Getting the full topology (segments level) from hoc

Anything that doesn't fit elsewhere.
hines
Site Admin
Posts: 1687
Joined: Wed May 18, 2005 3:32 pm

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

Post by hines »

The problem is that you are interating over the child sections when the number of children are changing within the body of the iterator.
(although it may not matter here and I can't tell if it does without thinking too hard about the code, it also does not seem to be a good idea to include
the Ext section within the forall processing). So I changed forall to 'forsec most' where most is the original sections. Anyway, here is the diff for the
changes I made.

Code: Select all

diff -r b66e93677626 test.hoc
--- a/test.hoc	Wed Nov 20 09:30:31 2013 -0500
+++ b/test.hoc	Wed Nov 20 09:40:31 2013 -0500
@@ -1,11 +1,20 @@
+    load_file("./sp.hoc")
+objref most
+most = new SectionList()
+forall most.append()
+
     objref somaref,secref,extVec,strFunc,stim,g1
     strFunc = new StringFunctions()
     create Exts[1]
-    proc hinesDisperseBranching(){localobj sl,secRef
+    proc hinesDisperseBranching(){local i localobj sl,secRef, clist
        counter = 0
-       forall{
+       forsec most{
           secRef = new SectionRef()
-          if(secRef.nchild>1){
+          if(secRef.nchild>2){
+  clist = new List()
+  for i=0, secRef.nchild-1 {
+    secRef.child[i] clist.append(new SectionRef())
+  }
              Exts[counter]{
                 L = 1e-9
                 diam = 1
@@ -18,7 +27,7 @@
                 for(x){
                    if (x > 0) {
                       printf("x is %f i is %d\n",x,i)
-                      connect secRef.child[i](0), x
+                      connect clist.o(i).sec(0), x
                       i +=1
                    }
                 }
@@ -34,14 +43,13 @@
     }
     proc countExts(){
     ext_num=0
-    forall{
+    forsec most {
        somaref= new SectionRef()
-       if (somaref.nchild>1){
+       if (somaref.nchild>2){
           ext_num+=1
        }
        }
     }
-    load_file("./sp.hoc")
     countExts()
     printf("Creating Exts %d\n",ext_num)
     create Exts[ext_num]

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

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

Post by roybens »

Thanks so much for this solution now the connection works well.
Next i set the NODEA and NODEB of the children using GETA and SETA in the branching.mod file as you recommended.
I am having difficulties now setting NODEA and NODEB of the Exts' nodes appropriately.
If i dont touch them they get a very large value since they have zero area and ressistance.
But i am not sure what value should i give them so they would result in 0 current.
Looking at triang and bksub in solve.c i thought that setting them to zero or to a very large number would result with NODERHS=0 but probably this is not the case.
Anyhow when i set them to zero fmatrix shows the updated values but after fcurrent the old values return. so i still cannot get the simulation to run as if there was no splitting.
I am attaching the relevant code
Thanks again for your time!
hoc file

Code: Select all

objref somaref,secref,extVec,strFunc,stim,g1,most
load_file("./sp.hoc")
most = new SectionList()
forall most.append()
strFunc = new StringFunctions()
create Exts[1]
proc hinesDisperseBranching(){local i localobj sl,secRef, clist, avec, bvec
	counter = 0
	forsec most{
		secRef = new SectionRef()
		if(secRef.nchild>2){
clist = new List()
  for i=0, secRef.nchild-1 {
    secRef.child[i] clist.append(new SectionRef())
  }
			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) {
						avec =new Vector()
						bvec =new Vector()
						printf("x is %f i is %d\n",x,i)
						clist.o(i).sec{
							avec =new Vector()
							bvec =new Vector()
							for(y){
								avec.append(GetA(y))
								bvec.append(GetB(y))
							}
						}
						connect clist.o(i).sec(0), x
						 k=0
					  clist.o(i).sec{
							for(y){
								SetA(y,avec.x(k))
								SetB(y,bvec.x(k))
							}
						}
						i +=1
					
					}
					
					
				}
				
			}
			fcurrent()
			secRef {
				connect Exts[counter](0), 1
			}
			Exts[counter]{
				for(x){ 
			printf("x is %f and sec is %s\n",x,secname())
					SetA(x,0)
					SetB(x,0)
					}
				}
			
	counter+=1		
		}
	}
}
proc countExts(){
ext_num=0
forsec most {
	somaref= new SectionRef()
	if (somaref.nchild>2){
		ext_num+=1
	}
	}
}

countExts()
printf("Creating Exts %d\n",ext_num)
fcurrent()
fmatrix()
create Exts[ext_num]
hinesDisperseBranching()
fmatrix()
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))
}
mod file

Code: Select all

COMMENT

ENDCOMMENT

NEURON { SUFFIX nothing }

VERBATIM
char* secname();
ENDVERBATIM
PROCEDURE init_files(){
	VERBATIM {
		
		
	}
	ENDVERBATIM
}





FUNCTION GetA(x) {
VERBATIM {
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif
Section* sec;
	Node* nd;
	sec = chk_access();
	if (_lx < 0. || _lx > 1.) {
	printf("_lx is %f and _lx*(double)(sec->nnode-1) is %f\n",_lx,_lx*(double)(sec->nnode-1));
		hoc_execerror("out of range, must be 0 < x <= 1", (char*)0);
	}
	if (_lx == 1.) {
		nd = sec->pnode[sec->nnode-1];
	}else{
		nd = sec->pnode[(int) (_lx*(double)(sec->nnode-1))];
	}
	return NODEA(nd);
}
ENDVERBATIM
}
FUNCTION GetB(x) {
VERBATIM {
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif
Section* sec;
	Node* nd;
	sec = chk_access();
	if (_lx < 0. || _lx > 1.) {
		printf("_lx is %f and _lx*(double)(sec->nnode-1) is %f\n",_lx,_lx*(double)(sec->nnode-1));
		hoc_execerror("out of range, must be 0 < x <= 1", (char*)0);
	}
	if (_lx == 1.) {
		nd = sec->pnode[sec->nnode-1];
	}else{
		nd = sec->pnode[(int) (_lx*(double)(sec->nnode-1))];
	}
	return NODEB(nd);
}
ENDVERBATIM
}
FUNCTION SetA(x,a) {
VERBATIM {
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif
Section* sec;
	Node* nd;
	sec = chk_access();
	if (_lx < 0. || _lx > 1.) {
		hoc_execerror("out of range, must be 0 < x <= 1", (char*)0);
	}
	if (_lx == 1.) {
		nd = sec->pnode[sec->nnode-1];
	}else{
		nd = sec->pnode[(int) (_lx*(double)(sec->nnode-1))];
	}
	printf("index is %d,NODEA(nd) is %f _la is %f\n",nd->v_node_index,NODEA(nd),_la);
	NODEA(nd) = _la;
}
ENDVERBATIM
}
FUNCTION SetB(x,b) {
VERBATIM {
#if defined(t)
	_NrnThread* _nt = nrn_threads;
#endif
Section* sec;
	Node* nd;
	sec = chk_access();
	if (_lx < 0. || _lx > 1.) {
		hoc_execerror("out of range, must be 0 < x <= 1", (char*)0);
	}
	if (_lx == 1.) {
		nd = sec->pnode[sec->nnode-1];
	}else{
		nd = sec->pnode[(int) (_lx*(double)(sec->nnode-1))];
	}
	printf("index is %d,NODEB(nd) is %f _lb is %f\n",nd->v_node_index,NODEB(nd),_lb);
	NODEB(nd) = _lb;
}
ENDVERBATIM
}

PROCEDURE MyPrintMatrix() {
VERBATIM {
	Section* sec;
	FILE* fm;
	fm= fopen("C:\fmatrix.dat", "wb");
	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
}
roybens
Posts: 54
Joined: Fri Mar 14, 2008 7:57 am

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

Post by roybens »

Hi again,
I am back in trying to apply this method on a more complicated model - one of the BBP pyramidal neurons https://bbp.epfl.ch/nmc-portal/downloads.
For some reason when applying this method to Mainen's model (1996) it works well but when applying to BBP it does not work at all and it is like everything is disconnected.
So in short this thread above describes how to modify neurons morphology to make sure there are no more than two childs on the same segment - this is done for efficiency when computing on GPU.
The approach is to add a very small section namely Exts where we need to change the branching. and connect it to the parent. Then children are connected to Exts and more segment are added to Exts to make sure every segment has only one child.
The properties of Exts are
L=1e-9
Diam = 1
cm = 1e-9
Ra = 1e-9
When i started troubleshooting why it does not work on BBP models i saw that i cannot check the parent:

Code: Select all

oc>access  cADpyr232_L5_TTPC1_0fb1ca4724[0].apic[118]
oc>objref sroy
oc>sroy = new SectionRef()
oc>sroy.parent
c:\nrn\bin\nrniv.exe: stack empty
 near line 280
 sroy.parent
            ^
Is this the accepted behavior?
There should be a parent according to topology()
Thanks in advance
Roy
hines
Site Admin
Posts: 1687
Joined: Wed May 18, 2005 3:32 pm

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

Post by hines »

You should check whether a section has a parent using sref.has_parent() before actually trying to reference the parent.
Also when many sections are connected to the same rootnode or parent section, it may be useful to use sref.has_trueparent() and sref.trueparent().
Parent and trueparent can differ when the parent's parent_connection is equal to the parent's section_orientation.
That is,

Code: Select all

oc>create a, b, c
oc>forall nseg=3
oc>connect b(0), a(1)
oc>connect c(0), b(0)
oc>topology()

|---|       a(0-1)
     `--|       b(0-1)
     `--|       c(0-1)

	1 
oc>objref sr
oc>c sr = new SectionRef()
oc>sr.parent print secname()
b
oc>sr.trueparent print secname()
a
Post Reply