collapse a dendritic tree
collapse a dendritic tree
Hoc code to collapse a dendritic tree
Collapses a dendritic tree into 3 compartments "proximal", "middle" and "distal".
http://cns.iaf.cnrs-gif.fr/files/COLLAPSE/
It collapses the dendritic tree into 3 compartments. I wonder if anyone really good with .hoc knows how to manipulate the code to instead collapse the tree into 20 compartments (or any other number of compartments specified by the user)? This is a great tool - but I think this manipulation could make it even better. I would be very grateful to anyone that thinks they can do it and takes the time to look.
--------------
bit of BACKGROUND to the algorithm used if you interested:
it does this: "collapsing the dendritic structure into a few compartments based on the conservation of axial resistance, a theme previously proposed by Bush and Sejnowski (1993).Bush–Sejnowski algorithm: Collapsing technique based on conserving Ra rather than the membrane surface area. This is done by making the cross-sectional area of the equivalent cylinder equal to the sum of the cross-sectional areas of all the dendrites represented by that equivalent cylinder. The length of the equivalent cylinder is just the average length of all the dendrites represented by the equivalent cylinder.
-----------------
Collapses a dendritic tree into 3 compartments "proximal", "middle" and "distal".
http://cns.iaf.cnrs-gif.fr/files/COLLAPSE/
It collapses the dendritic tree into 3 compartments. I wonder if anyone really good with .hoc knows how to manipulate the code to instead collapse the tree into 20 compartments (or any other number of compartments specified by the user)? This is a great tool - but I think this manipulation could make it even better. I would be very grateful to anyone that thinks they can do it and takes the time to look.
--------------
bit of BACKGROUND to the algorithm used if you interested:
it does this: "collapsing the dendritic structure into a few compartments based on the conservation of axial resistance, a theme previously proposed by Bush and Sejnowski (1993).Bush–Sejnowski algorithm: Collapsing technique based on conserving Ra rather than the membrane surface area. This is done by making the cross-sectional area of the equivalent cylinder equal to the sum of the cross-sectional areas of all the dendrites represented by that equivalent cylinder. The length of the equivalent cylinder is just the average length of all the dendrites represented by the equivalent cylinder.
-----------------
Still really struggling with this. So would really really appreciate some help. Below I have cut and paste 1) the original code and 2) my best stab at modifying it (which is not working code - there must be something wrong with it, although I cannot see what) - to collapse a dendritic tree into 4 compartments instead of 3. Although I actually want to collapse it into 40 compartments. Help!! I would be so, so grateful to anyone that can help as this is holding me back quite a lot.
1) THE ORIGINAL
(note that the orignial can also be found here:
http://cns.iaf.cnrs-gif.fr/files/COLLAPSE/
)
2) MY VERSION (not working)
1) THE ORIGINAL
(note that the orignial can also be found here:
http://cns.iaf.cnrs-gif.fr/files/COLLAPSE/
)
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF2 = 71 // define the cutoff between "middle" and "distal"
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(4,0)
cycdiam=new Vector(4,0)
cycnin=new Vector(4,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5){
if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
//Lbysurf_areabyasis BSC_Lbysurf_areabylong
Asurf=AmL*PI*Amdiam
Bsurf=BmL*PI*Bmdiam
AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
//Lbydiam_areabyasis Lbydiam_areabylong
AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
AmBmL=(AmL+BmL)/2
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
mrgL.set(secP, newL)
mrgri.set(secP, newri)
mrgdiam.set(secP, newdiam)
tosec.set(secA, secP)
tosec.set(secB, secP)
}
}
proc tagAcyc(){
tosec.set(secA, -$1)
tocyc.set(secA, $1)
}
proc tagBcyc(){
tosec.set(secB, -$1)
tocyc.set(secB, $1)
}
proc tagABcyc(){
tagAcyc($1)
tagBcyc($1)
}
proc cycmake(){local sec, extent
initto()
initcyc()
setmrg()
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {
tocyc.set(sec, -1)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF1 && secpathri.x(sec)< CUTOFF2) {
tocyc.set(sec, -2)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF1) {
tocyc.set(sec, -3)
}
}
for extent=1,3 {
getnextABP(extent)
shortnms()
while(secA>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,3 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,3 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,3 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,3 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,3 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,3 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,3 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
cycmake()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,3 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
2) MY VERSION (not working)
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF2 = 30 // define the cutoff between "middle" and "distal"
CUTOFF3 = 71
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(4,0)
cycdiam=new Vector(4,0)
cycnin=new Vector(4,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5){
if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
//Lbysurf_areabyasis BSC_Lbysurf_areabylong
Asurf=AmL*PI*Amdiam
Bsurf=BmL*PI*Bmdiam
AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
//Lbydiam_areabyasis Lbydiam_areabylong
AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
AmBmL=(AmL+BmL)/2
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
mrgL.set(secP, newL)
mrgri.set(secP, newri)
mrgdiam.set(secP, newdiam)
tosec.set(secA, secP)
tosec.set(secB, secP)
}
}
proc tagAcyc(){
tosec.set(secA, -$1)
tocyc.set(secA, $1)
}
proc tagBcyc(){
tosec.set(secB, -$1)
tocyc.set(secB, $1)
}
proc tagABcyc(){
tagAcyc($1)
tagBcyc($1)
}
proc cycmake(){local sec, extent
initto()
initcyc()
setmrg()
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF3 && secpathri.x(sec)< 1e10) {
tocyc.set(sec, -1)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {
tocyc.set(sec, -2)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF1 && secpathri.x(sec)< CUTOFF2) {
tocyc.set(sec, -3)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF1) {
tocyc.set(sec, -4)
}
}
for extent=1,4 {
getnextABP(extent)
shortnms()
while(secA>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,4 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,4 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,4 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,4 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,4 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,4 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,4 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
cycmake()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,4 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
-
- Posts: 2
- Joined: Tue May 23, 2006 8:39 pm
in cycmake(), this line looks suspicious:
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {
it probably needs to be something like:
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< CUTOFF3) {
==========
...... but i would suggest approaching your revisions in a new way:
1) make sure original code works okay for a variety of parameters
2) refactor the original code to use a Vector for CUTOFFs and use
the length of that vector when rewriting lines such as:
for extent=1,4 {
and in other places as called for.
3) in each small step of the refactoring process, make sure that, when n = 3, the modified code produces results that are identical to the original code.
4) after all refactoring is completed, debug for n != 3
http://en.wikipedia.org/wiki/Refactoring
best,
mike
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {
it probably needs to be something like:
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< CUTOFF3) {
==========
...... but i would suggest approaching your revisions in a new way:
1) make sure original code works okay for a variety of parameters
2) refactor the original code to use a Vector for CUTOFFs and use
the length of that vector when rewriting lines such as:
for extent=1,4 {
and in other places as called for.
3) in each small step of the refactoring process, make sure that, when n = 3, the modified code produces results that are identical to the original code.
4) after all refactoring is completed, debug for n != 3
http://en.wikipedia.org/wiki/Refactoring
best,
mike
thanks mike so much for your kind reply. I have changed the code to incorporate the line change you suggested. Unfortunately the code still doesn't work. dam!
I understand the alternative approach you suggested is superiour but I'm afraid I don't think I have the programming capability to do it. I am really only just hanging on by my finger tips doing these line changes as we are now. If I can get this code working for 4 as opposed to 3 - I think I will know how to just flesh it out manually (using cuts and pastes and little edits etc.) for 40 compartments. Then I will be very, very happy. Would you be kind enough to take a look at the newest, and still non-working, version of the code I have cut and paste below? I would be so, so grateful.
I understand the alternative approach you suggested is superiour but I'm afraid I don't think I have the programming capability to do it. I am really only just hanging on by my finger tips doing these line changes as we are now. If I can get this code working for 4 as opposed to 3 - I think I will know how to just flesh it out manually (using cuts and pastes and little edits etc.) for 40 compartments. Then I will be very, very happy. Would you be kind enough to take a look at the newest, and still non-working, version of the code I have cut and paste below? I would be so, so grateful.
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF2 = 30 // define the cutoff between "middle" and "distal"
CUTOFF3 = 71
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(4,0)
cycdiam=new Vector(4,0)
cycnin=new Vector(4,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5){
if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
//Lbysurf_areabyasis BSC_Lbysurf_areabylong
Asurf=AmL*PI*Amdiam
Bsurf=BmL*PI*Bmdiam
AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
//Lbydiam_areabyasis Lbydiam_areabylong
AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
AmBmL=(AmL+BmL)/2
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
mrgL.set(secP, newL)
mrgri.set(secP, newri)
mrgdiam.set(secP, newdiam)
tosec.set(secA, secP)
tosec.set(secB, secP)
}
}
proc tagAcyc(){
tosec.set(secA, -$1)
tocyc.set(secA, $1)
}
proc tagBcyc(){
tosec.set(secB, -$1)
tocyc.set(secB, $1)
}
proc tagABcyc(){
tagAcyc($1)
tagBcyc($1)
}
proc cycmake(){local sec, extent
initto()
initcyc()
setmrg()
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF3 && secpathri.x(sec)< 1e10) {
tocyc.set(sec, -1)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< CUTOFF3) {
tocyc.set(sec, -2)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=CUTOFF1 && secpathri.x(sec)< CUTOFF2) {
tocyc.set(sec, -3)
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF1) {
tocyc.set(sec, -4)
}
}
for extent=1,4 {
getnextABP(extent)
shortnms()
while(secA>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,4 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,4 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,4 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,4 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,4 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,4 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,4 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
cycmake()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,4 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
-
- Posts: 2
- Joined: Tue May 23, 2006 8:39 pm
Im not putting synapses onto my dendritic tree so this isn't so much of a worry for me. Am still really, really stuck with the coding. Haven't got past the non-working coding iteration in my last post. Any help would really make my weekend and would be so gratefully received.Mike Neubig wrote:the distributed nature of a dendritic arbor creates an uncoupled environment for synapses and membrane channels. collapsing a dendritic arbor into one with fewer branches, or a single tree, or a long cyclinder couples them.
Thanks so much to Dimitris who gave me this tip:
---------
Are you sure these vectors have the correct length?
proc initcyc(){
cycL=new Vector(4,0)
cycdiam=new Vector(4,0)
cycnin=new Vector(4,0)
}
---------
This got me to a version of the code working for 4 compartments as opposed to 3. I have cut and paste this working code below. HOWEVER, when I then tried to further modify the code to yield 10 compartments as opposed to 4 it didnt work. I have cut and paste this code below as well (in the next post). Stumped again. Dam! Would be so, so grateful for another pair of eyes. There must be a stupid mistake in there.
1) WORKING CODE FOR 4 COMPARTMENTS
---------
Are you sure these vectors have the correct length?
proc initcyc(){
cycL=new Vector(4,0)
cycdiam=new Vector(4,0)
cycnin=new Vector(4,0)
}
---------
This got me to a version of the code working for 4 compartments as opposed to 3. I have cut and paste this working code below. HOWEVER, when I then tried to further modify the code to yield 10 compartments as opposed to 4 it didnt work. I have cut and paste this code below as well (in the next post). Stumped again. Dam! Would be so, so grateful for another pair of eyes. There must be a stupid mistake in there.
1) WORKING CODE FOR 4 COMPARTMENTS
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF2 = 30 // define the cutoff between "middle" and "distal"
CUTOFF3 = 71
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(5,0)
cycdiam=new Vector(5,0)
cycnin=new Vector(5,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF3 && secpathri.x(sec)<1e10>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,4 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,4 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,4 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,4 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,4 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,4 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,4 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
cycmake()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,4 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
2) NON-WORKING CODE FOR 10 COMPARTMENTS
NON_WORKING - HELP!!!!!!
NON_WORKING - HELP!!!!!!
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF2 = 10 // define the cutoff between "middle" and "distal"
CUTOFF3 = 20
CUTOFF4 = 30
CUTOFF5 = 40
CUTOFF6 = 50
CUTOFF7 = 60
CUTOFF8 = 65
CUTOFF9 = 71
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(11,0)
cycdiam=new Vector(11,0)
cycnin=new Vector(11,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,10 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,10 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,10 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,10 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,10 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,10 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,10 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
cycmake()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
The error is upstream in line 358:mikey wrote:The error message returned is thus:
nrniv: syntax error in C:/address near line 651
I checked in the coding and this is line 651:
for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
Code: Select all
if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
} /* <- */
}
When deleting line 647 (a call to a non-existing proc/func cycmake), your code will be read without error messages.
Christoph
Crisoph, thanks so, so much for taking the time to look at this. I have made the 2 changes you suggested:
1) Take out the line:
cycmake()
2) Take out this bracket on line 358:
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
} I TOOK THIS BRACKET OUT!!!!
}
One thing I should mention about this bracket is that I didn't pick it up as an "orphan" bracket. It did have a corresponding opening bracket according to my Jedit text editor. But I took out this closing bracket regardless.
Upon making these 2 changes - the code did compile without any error messages and returned output. BUT the output was no good - because it gave the length and diameter of all the 10 compartments as 0. Which is obviously no good. So, Im still somewhat stuck. I have put the code with these changes below.
1) Take out the line:
cycmake()
2) Take out this bracket on line 358:
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
} I TOOK THIS BRACKET OUT!!!!
}
One thing I should mention about this bracket is that I didn't pick it up as an "orphan" bracket. It did have a corresponding opening bracket according to my Jedit text editor. But I took out this closing bracket regardless.
Upon making these 2 changes - the code did compile without any error messages and returned output. BUT the output was no good - because it gave the length and diameter of all the 10 compartments as 0. Which is obviously no good. So, Im still somewhat stuck. I have put the code with these changes below.
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF2 = 10 // define the cutoff between "middle" and "distal"
CUTOFF3 = 20
CUTOFF4 = 30
CUTOFF5 = 40
CUTOFF6 = 50
CUTOFF7 = 60
CUTOFF8 = 65
CUTOFF9 = 71
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(11,0)
cycdiam=new Vector(11,0)
cycnin=new Vector(11,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,10 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,10 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,10 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,10 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,10 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,10 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,10 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
Thanks so much to Dimitris who gave this tip:
"I have the impression that the problem is that you are trying to access memory which is not allocated, e.g. accessing beyond the bounds of an array. If you try a smaller number of compartments (e.g. 5 or 6) do you still get the problem?"
I have tried with 5 compartments - and I do still get this problem - not working.
"I have the impression that the problem is that you are trying to access memory which is not allocated, e.g. accessing beyond the bounds of an array. If you try a smaller number of compartments (e.g. 5 or 6) do you still get the problem?"
I have tried with 5 compartments - and I do still get this problem - not working.
I am so, so sorry. But there has been an awful mistake. The coding got corrupted in the cutting and pasting to the forum and the code for the working 4 compartments and the code for the non-working 10 compartments was NOT cut and paste correctly. SO, I now will cut and paste, and check, the code again to see that it gets replicated accurately in the forum. I am so, so sorry. I dont really know how it has happened.
UPDATE - just tried posting again, and the posting corrupts the code again. Will have to contact the administrators about this. Will hopefully get this sorted over next hour or so.
UPDATE - just tried posting again, and the posting corrupts the code again. Will have to contact the administrators about this. Will hopefully get this sorted over next hour or so.
When I wrote that the code would be read without errors after deleting the call to cycmake(), I didn't mean this wouldn't alter the behaviour of your code. Actually, cycmake() seems to be the central part of it.
Unfortunately, I had some problems reading your modified code, especially this if-clause:
Did you disable html in your post? Enabling html sometimes leads to strange behaviour when posting code. Some errors might have originated from copying & pasting html-spoiled code.
Anyway, I think it would be a great idea to use the excellent to-do-list Mike Neubig posted earlier. To get you going, I have started refactoring the code simply by introducing a CUTOFF vector instead of distinct variables. For readability puposes, I will only post the parts of the code that I have modified. Line numbers refer to the original code.
Line 72:
Line 415:
From there, I think it would be a better idea to use a variable such as n_comp instead of explicit numbers. For instance, instead of using
as you did you could use That way, you will only need to change the number of compartments once in your code. This can save you hours of debugging.
Christoph
Unfortunately, I had some problems reading your modified code, especially this if-clause:
Code: Select all
if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {
/*...*/
}
Anyway, I think it would be a great idea to use the excellent to-do-list Mike Neubig posted earlier. To get you going, I have started refactoring the code simply by introducing a CUTOFF vector instead of distinct variables. For readability puposes, I will only post the parts of the code that I have modified. Line numbers refer to the original code.
Line 72:
Code: Select all
// Modified 09/27/06, CSH------------------------------------------------
n_comp=3 //number of compartments
objectvar CUTOFF
CUTOFF=new Vector(n_comp-1)
CUTOFF.x[0] = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF.x[1] = 71 // define the cutoff between "middle" and "distal"
// End modified----------------------------------------------------------
Code: Select all
// Modified 09/27/06, CSH------------------------------------------------
for sec=1,NSEC-1 {
if(secpathri.x(sec)>=CUTOFF.x[n_comp-2] && secpathri.x(sec)< 1e10) {
tocyc.set(sec, -1)
}
}
for (i=n_comp-3;i>=0;i-=1) {
for sec=1,NSEC-1 {
if(secpathri.x(sec)>=CUTOFF.x[i] && secpathri.x(sec)< CUTOFF.x[i+1]) {
tocyc.set(sec, -n_comp+i+1)
}
}
}
for sec=1,NSEC-1{
if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF.x[0]) {
tocyc.set(sec, -n_comp)
}
}
// End modified ----------------------------------------------------------
Code: Select all
for extent=1,10
Code: Select all
for extent=1,n_comp
Christoph
Thanks so so much csh. Great!! What you have put looks great. What you have posted looks around about the same kind of thing that Michael Hines has just kindly sent me. What do you think? I think I'm right in saying your's and Michael's solutions are equivalent, but am not 100 percent. I have cut and paste the code that Michael Hines has sent me along with the explanatory email. The code should come out right now because I have turned HTML off now (excellent tip csh!!) - and from what I can see that solves the problem.
(NB. the Hines code below has been modified to yield 10 compartments. He originally sent me a version with for 3 compartments. Changed this line:
CUTOFF.append(0, 7.3, 71, 1e10)
to
CUTOFF.append(0, 7.3, 10, 20, 30, 40, 50, 60, 65, 71, 1e10)
---------
Enclosed is a collapse_n.hoc factorized as suggested by Mike Neubig
and based on the diff between the original collapse.hoc and your
putatively working 4 compartment variant. Collapse.hoc
and collapse_n.hoc produce identical results for 3 compartments.
I have no idea if the generalization to n compartments makes sense.
This is merely a code transformation that is identical to the
original for the CUTOFF.size == 4 case (three compartments)
where the first cutoff is 0 and the last cutoff is 1e10. If you wish to make additional compartments then just insert monotonically increasing numbers into
CUTOFF.append(0, 7.3, 71, 1e10)
You have my permission to add all my correspondence with regard to this
matter into the forum. In fact, please do so. Perhaps someone will
comment on the utility of the generalization. Also it is important
for the moderators to be aware of your experience with missing code
in the code fragments you posted.
------------
(NB. the Hines code below has been modified to yield 10 compartments. He originally sent me a version with for 3 compartments. Changed this line:
CUTOFF.append(0, 7.3, 71, 1e10)
to
CUTOFF.append(0, 7.3, 10, 20, 30, 40, 50, 60, 65, 71, 1e10)
---------
Enclosed is a collapse_n.hoc factorized as suggested by Mike Neubig
and based on the diff between the original collapse.hoc and your
putatively working 4 compartment variant. Collapse.hoc
and collapse_n.hoc produce identical results for 3 compartments.
I have no idea if the generalization to n compartments makes sense.
This is merely a code transformation that is identical to the
original for the CUTOFF.size == 4 case (three compartments)
where the first cutoff is 0 and the last cutoff is 1e10. If you wish to make additional compartments then just insert monotonically increasing numbers into
CUTOFF.append(0, 7.3, 71, 1e10)
You have my permission to add all my correspondence with regard to this
matter into the forum. In fact, please do so. Perhaps someone will
comment on the utility of the generalization. Also it is important
for the moderators to be aware of your experience with missing code
in the code fragments you posted.
------------
Code: Select all
/* ---------------------------------------------------------------------
HOC CODE TO COLLAPSE A DENDRITIC TREE
=====================================
Collapse a dendritic tree into 3 compartments "proximal", "middle"
and "distal".
The collapse is made such as the collapsed dendritic structure
preserves the axial resistance of the original structure.
The algorithm works by merging successive pairs of dendritic
branches into an equivalent branch (a branch that preserves the
axial resistance of the two original branches).
This merging of branches can be done according to different methods
selectable in the present code:
AREABYLONG
The code rejects the smallest original branch if it is more than
5 times smaller than the long orginal branch. If AREABYLONG is
set to 1, then the small branch is rescaled to the longest one
(a new diameter is assigned such that the small branch conserves
its axial resistance). This minimizes errors for dendritic
structures where branches have a large range of possible lengths.
AREABYLONG = 0 is the default.
AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF
During the merging of two branches, the length of the equivalent
branch is calculated by one out of three possible methods.
First, a simple average of the two branch length:
L = (L1+L2)/2
This is the method originally used by Bush and Sejnowski
(J. Neurosci. Methods 46: 159-166, 1993).
It is selected by setting AVGLENGTH=1.
Second, by weighting this average with the diameters of each branch
L = (diam1*L1 + diam2*L2)/(diam1+diam2)
This algorithm produces more fair merging in the case two branches of
very different diameters are to be merged. It is selected by setting
AVGLENGTHWDIAM=1.
Third, by weighting the average with membrane area:
L = (area1*L1 + area2*L2)/(area1+area2)
This algorithm is better in the case two branches of very different
areas are to be merged. It is selected by setting AVGLENGTHWSURF=1.
CUTOFF1, CUTOFF2
The dendritic branches are assigned to one of the three compartments
("proximal", "middle" or "distal") based on their path axial
resistance (in Meg-Ohms) to the soma. The cutoff values (CUTOFF1
and CUTOFF2) determine these cut-offs. For example, CUTOFF1=10 means
that all branches laying within 10 Meg-Ohm of axial resistance from
the soma will be merged together into the same ("proximal") section.
Written by M. Neubig & A. Destexhe
Laval University, 1997; CNRS 2000
-----------------------------------------------------------------------*/
objref CUTOFF
CUTOFF = new Vector()
CUTOFF.append(0, 7.3, 10, 20, 30, 40, 50, 60, 65, 71, 1e10)
n=CUTOFF.size
//CUTOFF1 = 7.3 // define the cutoff between "proximal" and "middle"
//CUTOFF2 = 71 // define the cutoff between "middle" and "distal"
AVGLENGTH = 0 // if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1 // if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0 // if 1, merge based on avg length, weighted by area
AREABYLONG = 0 // if 1, rescale the length before merging
xopen("tc-cell.geo") // read geometry file
Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0
Lbylong_areabylong = 0
if(AREABYLONG==1) {
print "Merge by rescaling length"
if(AVGLENGTH == 1) {
Lbyavg__areabylong = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabylong = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabylong = 1
}
} else {
print "Do not rescale length"
if(AVGLENGTH == 1) {
Lbyavg__areabyasis = 1
print "Merging method based on averaged length"
} else if(AVGLENGTHWDIAM == 1) {
Lbydiam_areabyasis = 1
print "Merging method based on averaged length, weigthed by diameters"
} else if(AVGLENGTHWSURF == 1) {
Lbysurf_areabyasis = 1
}
}
forall nseg=10
forall Ra=260
objectvar secaddress
objectvar paraddress
objectvar order
objectvar parsecndx
objectvar secri
objectvar secpathri
objectvar secL
objectvar secRa
objectvar secCm
objectvar mrgL
objectvar mrgdiam
objectvar mrgri
objectvar tosec
objectvar tocyc
objectvar rvp
objectvar slsoma
objectvar slroot
objectvar cycL
objectvar cycdiam
objectvar cycnin
proc initsec(){
secaddress=new Vector(NSEC,0) //set/reset
paraddress=new Vector(NSEC,0) //set/reset
order=new Vector(NSEC,0) //set/reset
parsecndx=new Vector(NSEC,0) //set/reset
secri=new Vector(NSEC,0) //set/reset
secpathri=new Vector(NSEC,0.0) //set/reset
secL=new Vector(NSEC,0) //set/reset
secRa=new Vector(NSEC,0) //set/reset
secCm=new Vector(NSEC,0) //set/reset
}
proc initmrg(){
mrgL=new Vector(NSEC,0) //set/reset
mrgdiam=new Vector(NSEC,0) //set/reset
mrgri=new Vector(NSEC,0) //set/reset
}
proc initcyc(){
cycL=new Vector(n,0)
cycdiam=new Vector(n,0)
cycnin=new Vector(n,0)
}
proc initto(){
tosec=new Vector(NSEC,-999) //set/reset
tocyc=new Vector(NSEC,-999) //set/reset
}
proc initVectors(){
initto()
initsec()
initmrg()
}
proc arborcharacterize(){local sec,sec1,sec2, aaa, x
{NSEC=0 forall NSEC=NSEC+1}
initVectors()
initcyc()
sec=0
forall {
if(sec==0){
secaddress.set(0, this_section())
paraddress.set(0, parent_section(0))
secRa.set(0,Ra)
secL.set(0,L)
secCm.set(0,cm)
mrgL.set(0,L)
{for(x) if(x>0) aaa=aaa+area(x) mrgdiam.set(0,aaa/PI/L)}
mrgri.set(0,secri.x(0))
}
if(sec!=0){
secaddress.set(sec, this_section())
paraddress.set(sec, parent_section(0))
{slsoma=new SectionList() rvp=new RangeVarPlot("v")
{soma rvp.begin(.5)} rvp.end(.5) rvp.list(slsoma)}
{slroot=new SectionList() rvp=new RangeVarPlot("v")
{ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)} ss=ss+1}}
rvp.end(.5) rvp.list(slroot)}
forsec slroot order.set(sec, order.x(sec)+1)
for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
if(secri.x(sec)>9999) secri.set(sec, -9999)
forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
}sec=sec+1}
setsec()
setmrg()
for sec1=1,NSEC-1 {
for sec2=0,NSEC-1 {
if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
}
}
}
proc setsec(){
sec=0
forall {
if(sec!=0){
secRa.set(sec,Ra)
secL.set(sec,L)
secCm.set(sec,cm)
}sec=sec+1}
}
proc setmrg(){
sec=0
forall {
if(sec!=0){
mrgL.set(sec,L)
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
mrgri.set(sec,secri.x(sec))
}sec=sec+1}
}
proc getnextABP(){local extent, sec
extent=$1
secA=0
for sec=1,NSEC-1 {if(tosec.x(sec)==-999){
if(order.x(sec)>=order.x(secA)){
if(abs(tocyc.x(sec))==extent){
secA=sec
}}}}
secB=0
for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){
if(sec!=secA){
if(tosec.x(sec)==-999){
if(abs(tocyc.x(sec))==extent){
secB=sec
}}}}}
secP=0
if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}
proc shortnms(){
{AmRa=0 AmL=0 Amdiam=0 Amri=0}
{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
if(secA>0) {
AmRa=secRa.x(secA)
AmL=mrgL.x(secA)
Amdiam=mrgdiam.x(secA)
Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
}
if(secB>0) {
BmRa=secRa.x(secB)
BmL=mrgL.x(secB)
Bmdiam=mrgdiam.x(secB)
Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
}
if(secP>0) {
PmRa=secRa.x(secP)
PmL=mrgL.x(secP)
Pmdiam=mrgdiam.x(secP)
Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
}
if(AREABYLONG==1){
if(secA>0 && secB>0){
if(AmL>BmL){
{BmL=AmL mrgL.set(secB,BmL)}
{Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
}
if(BmL>AmL){
{AmL=BmL mrgL.set(secA,AmL)}
{Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
}
}
}
}
proc joinAP(){local newRa, newL, newri, newdiam
newRa=(AmRa+PmRa)/2
newL=AmL+PmL
newri=Amri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
newRa=(BmRa+PmRa)/2
newL=BmL+PmL
newri=Bmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
mrgL.set(secP, newL)
mrgdiam.set(secP, newdiam)
mrgri.set(secP, newri)
tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam
AmL2BmL=AmL/BmL
BmL2AmL=BmL/AmL
if(AmL2BmL>=5) {
joinAP()
tosec.set(secB, -997)
}
if(BmL2AmL>=5) {
joinBP()
tosec.set(secA, -997)
}
if(AmL2BmL<5 && BmL2AmL<5){
if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
//Lbysurf_areabyasis BSC_Lbysurf_areabylong
Asurf=AmL*PI*Amdiam
Bsurf=BmL*PI*Bmdiam
AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
//Lbydiam_areabyasis Lbydiam_areabylong
AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
AmBmL=(AmL+BmL)/2
AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
AmBmcross=PI*(AmBmdiam^2)/4
AmBmRa=(AmRa+BmRa)/2
AmBmri=AmBmRa*AmBmL/AmBmcross/100
newRa=(AmBmRa+PmRa)/2
newL=AmBmL+PmL
newri=AmBmri+Pmri
newdiam=sqrt(newRa*newL*4/newri/PI/100)
}
mrgL.set(secP, newL)
mrgri.set(secP, newri)
mrgdiam.set(secP, newdiam)
tosec.set(secA, secP)
tosec.set(secB, secP)
}
}
proc tagAcyc(){
tosec.set(secA, -$1)
tocyc.set(secA, $1)
}
proc tagBcyc(){
tosec.set(secB, -$1)
tocyc.set(secB, $1)
}
proc tagABcyc(){
tagAcyc($1)
tagBcyc($1)
}
proc cycmake(){local ic, sec, extent
initto()
initcyc()
setmrg()
for ic = 2, n {
for sec = 1, NSEC-1 {
if (secpathri.x(sec) >= CUTOFF.x[n-ic] \
&& secpathri.x(sec) < CUTOFF.x[n - ic + 1]) {
tocyc.set(sec, 1 - ic)
}
}
}
for extent=1,n-1 {
getnextABP(extent)
shortnms()
while(secA>0) {
if (secB==0 && secP!=0) joinAP()
if (secB!=0 && secP!=0) joinABP()
if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
}
}
if(Lbylong_areabylong==1){
//Lbylong_areabylong
for extent=1,n-1 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabyasis==1){
//Lbysurf_areabyasis
for extent=1,n-1 {
totsurf=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbysurf_areabylong==1){
//Lbysurf_areabylong
for extent=1,n-1 {
totsurf=0
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
totsurf=totsurf+surf
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabyasis==1){
//Lbydiam_areabyasis
for extent=1,n-1 {
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbydiam_areabylong==1){
//Lbydiam_areabylong
for extent=1,n-1 {
maxL=0
totdiam=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
totdiam=totdiam+mrgdiam.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabyasis==1){
//Lbyavg__areabyasis
for extent=1,n-1 {
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
if(Lbyavg_areabylong==1){
//Lbyavg__areabylong
for extent=1,n-1 {
maxL=0
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
cycnin.set(extent, cycnin.x(extent)+1)
}
}
for sec=1,NSEC-1{
if(tocyc.x(sec)==extent) {
cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
}
}
if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
if(cycnin.x(extent)==0) cycL.set(extent, 0)
cycdiam.set(extent, sqrt(cycdiam.x(extent)))
}
}
}
arborcharacterize()
cycmake()
print "Section\tLength (um)\tDiameter (um)"
for ii=1,n-1 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)