Drosophila 3rd instar larval aCC motoneuron (Gunay et al. 2015)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:152028
Single compartmental, ball-and-stick models implemented in XPP and full morphological model in Neuron. Paper has been submitted and correlates anatomical properties with electrophysiological recordings from these hard-to-access neurons. For instance we make predictions about location of the spike initiation zone, channel distributions, and synaptic input parameters.
Reference:
1 . Günay C, Sieling FH, Dharmar L, Lin WH, Wolfram V, Marley R, Baines RA, Prinz AA (2015) Distal spike initiation zone location estimation by morphological simulation of ionic current filtering demonstrated in a novel model of an identified Drosophila motoneuron. PLoS Comput Biol 11:e1004189 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Drosophila;
Cell Type(s):
Channel(s): I Na,p; I Na,t; I A; I K;
Gap Junctions:
Receptor(s): Cholinergic Receptors;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; XPP; MATLAB;
Model Concept(s):
Implementer(s): Gunay, Cengiz [cgunay at emory.edu]; Sieling, Fred [fred.sieling at gmail.com]; Prinz, Astrid [astrid.prinz at emory.edu];
Search NeuronDB for information about:  Cholinergic Receptors; I Na,p; I Na,t; I A; I K;
/
Gunay_etal_2014
neuron-model
aCC-L3-neuron.hoc
aCC-L3-neuron+electrode.xml
aCC-L3-neuron-swc.hoc
calc-impedance.hoc
chan-DmKA-Marley.hoc
chan-DmKdr-Marley.hoc
chan-DmNaP-DmNav10.hoc
chan-DmNaT-ODowd.hoc
collapse-neuron-tree.hoc
current-inj-50pA-read-mV_dt_0.025ms.bin
data-axon-tail2-axon-50um-vc-noKdr-long-back-85mV-Na_4_lines_dt_0.025000ms.bin
data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-Na_4_lines_dt_0.025000ms.bin
data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-Na-5xNaP_4_lines_dt_0.025000ms.bin
data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-Na-5xNaT_4_lines_dt_0.025000ms.bin
data-axon-tail2-axon-70um-vc-noKdr-long-back-85mV-passive_4_lines_dt_0.025000ms.bin
data-axon-tail2-chans-axon_11_lines_dt_0.025000ms.bin
data-axon-tail2-chans-axon-last_11_lines_dt_0.025000ms.bin
data-axon-tail2-chans-botdend_11_lines_dt_0.025000ms.bin
data-axon-tail2-chans-ext-axon-70um_11_lines_dt_0.025000ms.bin
data-axon-tail2-chans-in-all_11_lines_dt_0.025000ms.bin
data-i-syn-10syns-20-EPSCs-10x-10ms-VC-60mV_6_lines_dt_0.025000ms.bin
data-i-syn-4dends-50-EPSCs-10x-10ms-VC-60mV_5_lines_dt_0.025000ms.bin
data-i-vclamp-syn-dend-513-180-EPSCs-10x-1ms-saturating_2_lines_dt_0.025000ms.bin
data-syn-dend-357_2_lines_dt_0.025000ms.bin
data-syn-dend-513_2_lines_dt_0.025000ms.bin
data-syn-dend-520_2_lines_dt_0.025000ms.bin
data-syn-dend-685_2_lines_dt_0.025000ms.bin
data-v-syn-10dends-20-EPSCs-10x-10ms-noVC_6_lines_dt_0.025000ms.bin
data-v-syn-4dends-50-EPSCs-10x-10ms-noVC_6_lines_dt_0.025000ms.bin
data-v-syn-dend-513-180-EPSCs-10x-1ms-saturating-noVC_5_lines_dt_0.025000ms.bin
data-v-syn-dend-685-AP_3_lines_dt_0.025000ms.bin
exp-axon-tail2.ses
exp-axon-tail2-chans-axon.ses
exp-axon-tail2-chans-axon-last.ses
exp-axon-tail2-chans-botdend.ses
exp-axon-tail2-chans-ext-axon-50um-onlyNa.ses
exp-axon-tail2-chans-ext-axon-70um.ses
exp-axon-tail2-chans-ext-axon-70um-10alphasynapses.ses
exp-axon-tail2-chans-ext-axon-70um-10x-mimic-sustained.ses
exp-axon-tail2-chans-ext-axon-70um-10x-mimic-sustained-random.ses
exp-axon-tail2-chans-ext-axon-70um-mimic-synapses.ses
exp-axon-tail2-chans-ext-axon-70um-mimic-synapses-sustained-currents.ses
exp-axon-tail2-chans-ext-axon-70um-mimic-synapses-v-change.ses
exp-axon-tail2-chans-ext-axon-70um-onlyNa.ses
exp-axon-tail2-chans-ext-axon-70um-tomasz.ses
exp-axon-tail2-chans-in-all.ses
figures.m
fitfuncs.hoc
graph-i-vc-ext-axon.ses
iclamp-50pA.ses
IClamp-steps.ses
inc-first.ses
lincir-vclamp.hoc
lincir-vclamp.ses
NaP_NaT_data.csv
neuron-CB.ses
neuron-CB+electrode.hoc
neuron-CB-act-electrode-embed-IClamp.ses
neuron-CB-ext-axon.ses
neuron-CB-ext-axon-2pieces.ses
neuron-CB-ext-axon-2pieces-chans-axon.ses
neuron-CB-ext-axon-2pieces-chans-axon-last.ses
neuron-CB-ext-axon-2pieces-chans-botdend.ses
neuron-CB-ext-axon-2pieces-chans-ext-axon-50um-onlyNa.ses
neuron-CB-ext-axon-2pieces-chans-ext-axon-70um.ses
neuron-CB-ext-axon-2pieces-chans-ext-axon-70um-10alphasynapses.ses *
neuron-CB-ext-axon-2pieces-chans-ext-axon-70um-10x-mimic-sustained.ses *
neuron-CB-ext-axon-2pieces-chans-ext-axon-70um-mimic-synapses.ses *
neuron-CB-ext-axon-2pieces-chans-ext-axon-70um-mimic-synapses-v-change.ses *
neuron-CB-ext-axon-2pieces-chans-ext-axon-70um-onlyNa.ses
neuron-CB-ext-axon-2pieces-chans-in-all.ses
neuron-CB-pas-electrode-embed.ses
neuron-CB-pas-electrode-embed-fit-pas.ses
neuron-CB-pas-electrode-embed-fit-pas-VClamp.ses
neuron-CB-pas-electrode-embed-IClamp.ses
neuron-CB-pas-electrode-embed-test-axon-hh-chans.ses
neuron-Import3D-CellBuilder.ses
neuron-NL-CellBuilder.ses
neuron-NL-CellBuilder-pas.ses
neuron-NL-CellBuilder-pas-electrode.ses
neuron-NL-CellBuilder-pas-Na.ses
neuron-PointProcessMgr-ext-axon-2pieces-chans-ext-axon-70um-10alphasynapses.ses
nrn-fit-cap-02_dt_0.025000ms_dy_1e-9nA.bin
shape-plot.ses
SkeletonTree_ORR_aCC_48h1_NL.hoc
soma-vclamp-testbed.ses
stats.hoc
vclamp_-85_to_-25mV.ses
vclamp_soma_-60mV.ses
vclamp_soma_-60mV_syn1234.ses
vclamp_soma_-60mV_syni.ses
vclamp-family.ses
v-graph.ses
v-graph-bigger.ses
v-graph-bigger-axon-2pieces.ses
                            
/* ---------------------------------------------------------------------

	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)