Ca1 pyramidal neuron: reduction model (Marasco et al. 2012)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:146376
"... Here we introduce a new, automatic and fast method to map realistic neurons into equivalent reduced models running up to >40 times faster while maintaining a very high accuracy of the membrane potential dynamics during synaptic inputs, and a direct link with experimental observables. The mapping of arbitrary sets of synaptic inputs, without additional fine tuning, would also allow the convenient and efficient implementation of a new generation of large-scale simulations of brain regions reproducing the biological variability observed in real neurons, with unprecedented advances to understand higher brain functions."
Reference:
1 . Marasco A, Limongiello A, Migliore M (2012) Fast and accurate low-dimensional reduction of biophysically detailed neuron models. Sci Rep 2:928 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Dendrite;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampus CA1 pyramidal GLU cell;
Channel(s): I Na,t; I A; I K; I h;
Gap Junctions:
Receptor(s): AMPA;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Simplified Models; Detailed Neuronal Models;
Implementer(s): Limongiello, Alessandro [alessandro.limongiello at unina.it];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal GLU cell; AMPA; I Na,t; I A; I K; I h;
/
reduction1.0
MaxStim_Output
morphologies
Readme.html
distr.mod *
h.mod *
kadist.mod *
kaprox.mod *
kdrca1.mod *
na3n.mod *
naxn.mod *
clusterisingMethods.hoc
fixnseg.hoc *
MaxStimPROCEDURE1.0.hoc
mergingMethods.hoc
mosinit.hoc
ranstream.hoc *
REDUCTION1.0.hoc
screenshot.png
SoftReduction1.0.doc
stimulation1.hoc
useful&InitProc.hoc
                            
/*---------------------------------------------
                   USEFUL PROCEDURES
  ---------------------------------------------*/

func max(){
    if ($1 > $2){ 
       return $1
    } else {
      return $2
    }
}


objectvar factGmaxSynMorphVect

factGmaxSynMorphVect = new Vector(100,1.0)		//set/reset

factGmaxSynMorphVect.x(0)= 340
factGmaxSynMorphVect.x(1)= 180
factGmaxSynMorphVect.x(2)= 160
factGmaxSynMorphVect.x(3)= 220
factGmaxSynMorphVect.x(4)= 200
factGmaxSynMorphVect.x(5)= 200
factGmaxSynMorphVect.x(6)= 200
factGmaxSynMorphVect.x(7)= 200
factGmaxSynMorphVect.x(8)= 240
factGmaxSynMorphVect.x(9)= 240
factGmaxSynMorphVect.x(10)= 180  
factGmaxSynMorphVect.x(11)= 160  
factGmaxSynMorphVect.x(12)= 200  
factGmaxSynMorphVect.x(13)= 180  
factGmaxSynMorphVect.x(14)= 200  
factGmaxSynMorphVect.x(15)= 220  
factGmaxSynMorphVect.x(16)= 200 

strdef str,strFile

proc checkFile(){
     if (strcmp(strFile,"./morphologies/geoc70963.hoc")==0) {iMorph = 0
     }else if (strcmp(strFile,"./morphologies/geoc81462.hoc")==0) {iMorph = 1
     }else if (strcmp(strFile,"./morphologies/geoc70863.hoc")==0) {iMorph = 2
     }else if (strcmp(strFile,"./morphologies/geoc20466.hoc")==0) {iMorph = 3
     }else if (strcmp(strFile,"./morphologies/geocd2351.hoc")==0) {iMorph = 4
     }else if (strcmp(strFile,"./morphologies/geoc8076.hoc")==0) {iMorph = 5   
     }else if (strcmp(strFile,"./morphologies/geoc9236.hoc")==0) {iMorph = 6
     }else if (strcmp(strFile,"./morphologies/geoc91665.hoc")==0) {iMorph = 7
     }else if (strcmp(strFile,"./morphologies/geoc80761.hoc")==0) {iMorph = 8
     }else if (strcmp(strFile,"./morphologies/geocd0351.hoc")==0) {iMorph = 9 
     }else if (strcmp(strFile,"./morphologies/geoc30465.hoc")==0) {iMorph = 10    
     }else if (strcmp(strFile,"./morphologies/geocd1152.hoc")==0) {iMorph = 11
     }else if (strcmp(strFile,"./morphologies/geoc73162.hoc")==0) {iMorph = 12
     }else if (strcmp(strFile,"./morphologies/geoc72965.hoc")==0) {iMorph = 13
     }else if (strcmp(strFile,"./morphologies/geoc91662.hoc")==0) {iMorph = 14
     }else if (strcmp(strFile,"./morphologies/geoc62564mod.hoc")==0) {iMorph = 15
     }else if (strcmp(strFile,"./morphologies/geoc73166.hoc")==0) {iMorph = 16
     }else {iMorph = 99}
    
     maxStim = factGmaxSynMorphVect.x(iMorph)
     if (iMorph == 99) maxStimOnOff = 0
}


objref f

proc readingFile(){

    f = new File()
    f.chooser("", "Reading Morphologic Model", "", "Accept", "", ".\\morphologies")
    if (f.chooser()) {
       firstInit = 0
       f.getname(strFile)
       checkFile()
       
    }else firstInit = 1
    
    

}

proc AleFunctionalSetting(){
     
     if (optPrint==1) print "AleFunctionalSetting Setting"
     forall delete_section() 
     xopen(strFile)		// read geometry file
     firstInit = 0

     
     LENGTHMERGINGserialMETHOD = -1
     LENGTHMERGINGparallelMETHOD = 2
     DIAMMERGINGserialMETHOD = 0
     DIAMMERGINGparallelMETHOD = 0
     RESCALEMOD = 0
     CUTTINGFACT = 5
     PRESERVINGsomaMETHOD = 0
     PRESERVINGMETHOD = 0 
     RaMERGINGMETHOD=1
     SURFFACT_PARTIAL_TOTAL = 1
     SURFFACT_PARTIAL_METHOD=2
     withLocation = 1
     
     Ra0 = 150
     RaAx = 50
     g_pas0 = 1/28000                //Siemens/cm^2   1/30000 
     cm0 = 1

     forall {insert pas e_pas=Vrest g_pas=g_pas0 Ra=Ra0 cm=cm0}
     forsec "axon[0]" {insert pas e_pas=Vrest g_pas = g_pas0 Ra=RaAx cm=cm0}  

     geom_nseg()
     totOr=0
     forall {totOr=totOr+nseg}   
     if(optPrint) print "totSegIni: ", totOr 
}



/*---------------------------------------------
                   VARIABLE DEFINITION
  ---------------------------------------------*/
	
objectvar order			// number of sections between each section and soma 
objectvar parsec        // parentsection id for each section

// section data input
objectvar secri			
objectvar secpathri		
objectvar secpathri0
objectvar secpathImp		
objectvar secpathImp0
objectvar secL			
objectvar secRa			
objectvar secCm	
objectvar secdiam
objectvar secSurf
objectvar secApicalBasal
// section data during the merging		
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgdiamSer
objectvar mrgdiamPar
objectvar mrgri	
objectvar mrgri2
objectvar mrgSurf
// equivalent section data after the merging		
objectvar eqL		
objectvar eqdiam	
objectvar eqdiamSer
objectvar eqdiamPar
objectvar eqRa


objectvar surfFactVect
objectvar eqri2
objectvar eqsecpathri
objectvar eqsecpathri0
objectvar eqsecpathImp
objectvar eqsecpathImp0
objectvar orMinpathri
objectvar orMaxpathri
//objectvar eqri	
objectvar eqSurf
objectvar orSurf	
objectvar orSurf1
objectvar orSurf2
objectvar surfFact
objectvar surfFactRmCm

// cluster data
objectvar clusterParent    // link between each cluster and its parent
objectvar clusterParentPos // position of each cluster as regards its parent
objectvar clusterList      // List of clusters
objectvar newClusterList   // List of clusters after merging Y_group_sections
objectvar clusterMarker    // Marker to identify each section belonging to its cluster
objectvar somaVector       // Vector of sections belonging to the soma
objectvar axonVector       // Vector of sections belonging to the axon


	
	
/*---------------------------------------------
                   INIT PROCEDURE
  ---------------------------------------------*/	
	
	
proc initDataStruct(){

	order = new Vector(NSEC,0)		//set/reset
	parsec = new Vector(NSEC,0)		//set/reset

	secri = new Vector(NSEC,0)		//set/reset
	secpathri = new Vector(NSEC,0.0)		//set/reset
	secpathri0 = new Vector(NSEC,0.0)		//set/reset
	secpathImp = new Vector(NSEC,0.0)		//set/reset
	secpathImp0 = 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
	secdiam = new Vector(NSEC,0)		//set/reset
	secSurf = new Vector(NSEC,0)		//set/reset
	secApicalBasal = new Vector(NSEC,0)		//set/reset if 1=Apical elseif 0=Basal
	
	mrgL = new Vector(NSEC,0)   		//set/reset
	mrgdiam = new Vector(NSEC,0)		//set/reset
	mrgdiamSer = new Vector(NSEC,0)		//set/reset
	mrgdiamPar = new Vector(NSEC,0)		//set/reset	
	mrgri = new Vector(NSEC,0)		//set/reset
	mrgri2 = new Vector(NSEC,0)		//set/reset
	mrgSurf = new Vector(NSEC,0)		//set/reset
	
    clusterList = new List()	
	newClusterList = new List()
	clusterMarker = new Vector(NSEC,0)		//set/reset
	somaVector = new Vector(0,0)	
    axonVector = new Vector(0,0)	

}


/*---------------------------------------------
                   CELL ANALYSIS
  ---------------------------------------------*/	

// Useful variables
objectvar secNameList
objectvar secName
secNameList=new List()
objectvar secaddress	// section address	
objectvar paraddress	// parentsection address
// data to compute path from each section to the soma
objectvar rvp			
objectvar slsoma		
objectvar slroot


proc cellAnalysis(){local sec,sec1,sec2, x

	strdef sname

	{NSEC=0  forall if (!issection("eqSection.*")){NSEC=NSEC+1}}		
    initDataStruct()
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
            
	sec = 0
	forall if (!issection("eqSection.*")){
        sname =secname() secName = new String(sname) secNameList.append(secName)
        secaddress.set(sec, this_section())
        paraddress.set(sec, parent_section(0))
        secRa.set(sec,Ra)
        secL.set(sec,L)
        secCm.set(sec,cm) 
        secdiam.set(sec,diam)   
        mrgL.set(sec,L)
        if (y3d(1)>=2) secApicalBasal.set(sec,1)
        if (issection("soma.*")){
              for(x) if(x>0) {
                     secSurf.set(sec,secSurf.x(sec)+area(x))
                     secri.set(sec, secri.x(sec)+ri(x))
              }
              //mrgdiam.set(sec,aaa/PI/L)
              //mrgri.set(sec,secri.x(0))
              //print sec,totSomaArea,mrgL.x(sec),mrgdiam.x(sec),secri.x(sec)
              somaVector.append(sec)
        } else if (issection("axon[0]")){
              if(optPrint) print "pippo: ",sec
              for(x) if(x>0) {
                     secSurf.set(sec,secSurf.x(sec)+area(x))
                     secri.set(sec, secri.x(sec)+ri(x))
              }
              axonVector.append(sec)
        } else {
    		  {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) /*print sec,order.x(sec) psection()*/}
              for(x) if(x>0) {
                     secSurf.set(sec,secSurf.x(sec)+area(x))
                     secri.set(sec, secri.x(sec)+ri(x))
                     /*print sec,secri.x(sec)*/
              }
    		  if(secri.x(sec)>9999) secri.set(sec, -9999)
    		  nsecPath = 0
              forsec slroot {   nsecPath=nsecPath+1 
                                for(x) if(x>0){secpathri.set(sec,secpathri.x(sec)+ri(x))}
                            }
              nsecPath1 = 0
              forsec slroot {
                     nsecPath1=nsecPath1+1 
                     for(x) if(x>0 && nsecPath1<nsecPath) secpathri0.set(sec,secpathri0.x(sec)+ri(x))
              }

              if (PRESERVINGMETHOD==0){             //PRESERVINGMETHOD = 0 by Destexhe
                     mrgdiamSer.set(sec,secdiam.x(sec))
                     mrgdiamPar.set(sec,secdiam.x(sec)^2)
                     //print mrgdiamPar.x(sec)
                     mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))                     
              }else if (PRESERVINGMETHOD==1){
                     mrgdiam.set(sec,secSurf.x(sec)/secL.x(sec)/PI)
              }
              mrgri.set(sec,secri.x(sec))		  
              mrgri2.set(sec,secri.x(sec))	
              mrgSurf.set(sec, secSurf.x(sec))
        }
    sec=sec+1}	
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsec.set(sec1,sec2)}
		}
	}
}

/*---------------------------------------------
                   CLUSTERISING PROCEDURES
  ---------------------------------------------*/
xopen("clusterisingMethods.hoc")

/*---------------------------------------------
                   MERGING PROCEDURES
  ---------------------------------------------*/
xopen("mergingMethods.hoc")

/*---------------------------------------------
                   REDUCTION PROCEDURES
  ---------------------------------------------*/


// Useful variables
objectvar searchAvailable
objectvar mergeAvailable



// seaching the next Y_group_sections / output:  {secA,secB,secP} 
proc getNextY(){local ii, nsect, idA, idB, sec   //input: Vector of id sections

    // seaching the last available section
	secA=0			
	nsect = $o1.size()-1
	idA = 0 
	for ii=0,nsect {
        if(searchAvailable.x(ii)){
        	sec = $o1.x(ii)	
            if(order.x(sec)>=order.x(secA)) {secA=sec idA = ii}
        }
    }
    
    // seaching the brother section of secA
	secB=0
	idB = 0
	for ii=0,nsect {
        if(searchAvailable.x(ii)){
        	sec = $o1.x(ii)		
            if(parsec.x(sec)==parsec.x(secA)&&sec!=secA) {secB=sec idB = ii}
        }
    }

    // seaching the parent section of secA
	secP=0
	if ($o1.contains(parsec.x(secA))) secP=parsec.x(secA)	

    // marking unavailable sections already visited	
	searchAvailable.set(idA,0)
	if (secB!=0) searchAvailable.set(idB,0)
	
    if (secA!=0){
    	if (secP==0){ 
           mergeAvailable.set(idA,1)
    	   if (secB!=0) mergeAvailable.set(idB,1)
        }
    }
}


// Useful proc used in reduction() after each call of getNextY() and before each call of merge2() or mergeY()
//        input:   {secA,secB,secP}
//        output:  {AmRa,AmL,Amdiam,Amri,BmRa,BmL,Bmdiam,Bmri,PmRa,PmL,Pmdiam,Pmri}
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0 Amsurf=0 AmdiamSer=0 AmdiamPar=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0 Bmsurf=0 BmdiamSer=0 BmdiamPar=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0 Pmsurf=0 PmdiamSer=0 PmdiamPar=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA)
		AmdiamSer=mrgdiamSer.x(secA)
        AmdiamPar=mrgdiamPar.x(secA)
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
		Amsurf = AmL*Amdiam*PI
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		BmdiamSer=mrgdiamSer.x(secB)
        BmdiamPar=mrgdiamPar.x(secB)
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
		Bmsurf = BmL*Bmdiam*PI
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		PmdiamSer=mrgdiamSer.x(secP)
        PmdiamPar=mrgdiamPar.x(secP)
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
		Pmsurf = PmL*Pmdiam*PI
	}
    
	if(RESCALEMOD==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
          if (PRESERVINGMETHOD==0){             //PRESERVINGMETHOD = 0 by Destexhe
               Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100))
        	   BmdiamSer=Bmdiam
               BmdiamPar=Bmdiam
          }else if (PRESERVINGMETHOD==1){
               Bmdiam=Bmsurf/BmL/PI } 
          mrgdiam.set(secB,Bmdiam) 
          {Bmsurf = BmL*Bmdiam*PI}
	    }else if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      if (PRESERVINGMETHOD==0){             //PRESERVINGMETHOD = 0 by Destexhe
               Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100))
        	   AmdiamSer=Amdiam
               AmdiamPar=Amdiam               
          }else if (PRESERVINGMETHOD==1){
               Amdiam=Amsurf/AmL/PI } 
          mrgdiam.set(secA,Amdiam)
	      {Amsurf = AmL*Amdiam*PI}
	    }
	  }
	}
}


// merging procedure of the section secA with its parent secP
proc mergingParentSec(){
    mergingSerialMethod(AmRa,PmRa,AmL,PmL,Amdiam,Pmdiam,Amri,Pmri,Amsurf,Pmsurf,AmdiamSer,PmdiamSer,AmdiamPar,PmdiamPar)
    mrgL.set(secP, newL)
	mrgri.set(secP, newri)
	mrgri2.set(secP, Amri+Pmri)
	mrgSurf.set(secP, max(Amsurf,Pmsurf))
	mrgdiam.set(secP, newdiam) 
	mrgdiamSer.set(secP, newdiamSer) 
	mrgdiamPar.set(secP, newdiamPar) 
	//print "newdiamPar",newdiamPar,"secP",secNameList.o(secP).s
}


// merging procedure of Y_group_sections
proc mergingYSec(){local AmL2BmL,BmL2AmL
	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=CUTTINGFACT){mergingSerialMethod(AmRa,PmRa,AmL,PmL,Amdiam,Pmdiam,Amri,Pmri,Amsurf,Pmsurf,AmdiamSer,PmdiamSer,AmdiamPar,PmdiamPar)}
    if(BmL2AmL>=CUTTINGFACT){mergingSerialMethod(BmRa,PmRa,BmL,PmL,Bmdiam,Pmdiam,Bmri,Pmri,Bmsurf,Pmsurf,BmdiamSer,PmdiamSer,BmdiamPar,PmdiamPar)} 
	if(AmL2BmL<CUTTINGFACT && BmL2AmL<CUTTINGFACT){
       mergingYMethod(AmRa,BmRa,PmRa,AmL,BmL,PmL,Amdiam,Bmdiam,Pmdiam,Amri,Bmri,Pmri,AmdiamSer,BmdiamSer,PmdiamSer,AmdiamPar,BmdiamPar,PmdiamPar)}
	
    mrgL.set(secP, newL)
	mrgri.set(secP, newri)
	mrgri2.set(secP, newri2)
	mrgSurf.set(secP, newSurf)
	mrgdiam.set(secP, newdiam)  
	mrgdiamSer.set(secP, newdiamSer) 
	mrgdiamPar.set(secP, newdiamPar) 	
    //print "newdiamPar da Y",newdiamPar,"secP",secNameList.o(secP).s	
}

// Useful variables
objectvar clusterVectApp
objectvar mergeAvailableApp

// procedure to create a new Cluster List after merging Y_group_sections
proc newClusters(){local nsect, sec
    {clusterVectApp = $o1 mergeAvailableApp = $o2}
    nsect = clusterVectApp.size()
    clusterVect = new Vector(0,0)

	for ii=0,nsect-1 {
        if(mergeAvailableApp.x(ii)) {/*print clusterVectApp.x(ii)*/ clusterVect.append(clusterVectApp.x(ii))}
    }    
    newClusterList.append(clusterVect)

}

objectvar appVect

proc reduction(){local i, ii, nsect

	eqL=new Vector(clusterList.count,0)   		//set/reset
	eqdiam=new Vector(clusterList.count,0)		//set/reset
	eqdiamSer=new Vector(clusterList.count,0)		//set/reset
	eqdiamPar=new Vector(clusterList.count,0)		//set/reset	
	//eqRa=new Vector(clusterList.count,0)		//set/reset	
	eqri2=new Vector(clusterList.count,0)		//set/reset	
	eqSurf=new Vector(clusterList.count,0)		//set/reset	
	orSurf=new Vector(clusterList.count,0)		//set/reset
	orSurf1=new Vector(clusterList.count,0)		//set/reset
	orSurf2=new Vector(clusterList.count,0)		//set/reset
	surfFact=new Vector(clusterList.count,1)		//set/reset
	surfFactRmCm=new Vector(clusterList.count,1)		//set/reset	
	//eqri=new Vector(clusterList.count,0)		//set/reset
	
	
	newClusterList.append(somaVector)
	mergingSoma(newClusterList.o(0))
	iniClusterList = 1
	if (axonVector.size()>0){
    	newClusterList.append(axonVector)
    	mergingAxon(newClusterList.o(1))
    	iniClusterList = 2
    }
    	
	for i=iniClusterList,clusterList.count-1 {
        nsect = clusterList.o(i).size()   
        searchAvailable = new Vector(nsect,1)
        mergeAvailable = new Vector(nsect,0)
		getNextY(clusterList.o(i))
		shortnms()
		while(secA>0) {	
			if (secB==0 && secP!=0) mergingParentSec()
			if (secB!=0 && secP!=0) mergingYSec()
            //print secNameList.o(secA).s,secA,secNameList.o(secB).s,secB,secNameList.o(secP).s,secP
			getNextY(clusterList.o(i))
			shortnms()
		}
		newClusters(clusterList.o(i),mergeAvailable)
		mergingAllMethod(newClusterList.o(i),i,nsect)

        totSurf = 0
        appVect = new Vector(0,0)
        for ii=0,nsect-1 {
            appVect.append(secSurf.x(clusterList.o(i).x(ii)))
            totSurf = totSurf + secSurf.x(clusterList.o(i).x(ii))
        }
        //print totSurf
        
        if (SURFFACT_PARTIAL_METHOD==0 || SURFFACT_PARTIAL_METHOD==3){
            orSurf.set(i,appVect.sum())
            surfFact.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
            surfFactRmCm.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
        }else if (SURFFACT_PARTIAL_METHOD==1 || SURFFACT_PARTIAL_METHOD==4){
            orSurf.set(i,orSurf1.x(i))
            surfFact.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))
            surfFactRmCm.set(i,surfFactVect.x(i)*orSurf.x(i)/eqSurf.x(i))    
            //print surfFact.x(i),surfFactRmCm.x(i)
        }        
  
              
        		
	}
}



/*---------------------------------------------
            OUTPUT DATA:
            PRINT CLUSTERLIST
            CREATING EQUIVALENT CELL
  ---------------------------------------------*/

// Useful variables
objectvar clusterListApp

proc printClusterList(){local nCluster,secApp,i,j,ii
    clusterListApp = $o1
    nCluster = clusterListApp.count
    
    for i=0,nCluster-1{
        print "_____________________"
        print "cluster n. ", i
        
        for j=0,clusterListApp.o(i).size()-1{
            secApp = clusterListApp.o(i).x(j)
            
            print secApp, secNameList.o(secApp).s, secL.x(secApp),secdiam.x(secApp),mrgdiamPar.x(secApp), secSurf.x(secApp),secpathri0.x(secApp),secpathri.x(secApp)
        }
        
    }    
    print "Section\tLength (um)\tDiameter (um)\tOld Area (um^2)\tNew Area (um^2)\tArea Fact"
    
    for ii=0,nCluster-1 print ii, "\t", eqL.x(ii),"\t", eqdiam.x(ii),"\t",orSurf.x(ii),"\t",eqSurf.x(ii),"\t",surfFactRmCm.x(ii)

    print "eqRa\t orMinpathri\t orMaxpathri\t eqri2\t eqriEffettiva\t eqpathri0\t eqpathri"
    
    for ii=0,nCluster-1 {
        eqSection[ii] {
             if(RaMERGINGMETHOD==0){RaEq = factEqRa*eqRa.x(ii)}\
             else if(RaMERGINGMETHOD==1){RaEq = factEqRa*eqri2.x(ii)*PI*(diam/2)^2*100/L}\
             else if(RaMERGINGMETHOD==2){RaEq = factEqRa*(orMaxpathri.x(ii)-orMinpathri.x(ii))*PI*(diam/2)^2*100/L}
        }
        print RaEq,"\t",orMinpathri.x(ii),"\t",orMaxpathri.x(ii),"\t",eqri2.x(ii),"\t",(eqRa.x(ii)*eqL.x(ii))/(PI*((eqdiam.x(ii)^2)/4)*100),"\t",eqsecpathri0.x(ii),"\t",eqsecpathri.x(ii)
        }
    
     
}

// Useful variables
objectvar sh
objectvar shPlot
objref vbox0, vbox1

strdef app
objref strobj,rr,gg,bb
strobj = new StringFunctions()

proc coloringClusterList(){local nCluster,sec,col,i
    clusterListApp = $o1
    nCluster = clusterListApp.count

    sec = 0
    forall if (!issection("eqSection.*")){
           col = 1
           for (i = 0; i < nCluster; i = i + 1){
               //print "pippo",i,col,col-int(col/10)
               if (col-int(col/10)*10==2 || col-int(col/10)*10==0) col = col-int(col/10)*10+1
               //print "pluto",i,col
               if (clusterListApp.o(i).contains(sec)&&!issection(app)) sh.color(col)
               col = col + 1
           }
           sec=sec+1
    }
    
    col = 1
    forsec "eqSection.*" {if (col-int(col/10)*10==2 || col-int(col/10)*10==0) {col = col-int(col/10)*10+1} sh.color(col) col=col+1}    
    
    length = strobj.len(app)
    if (length>0) forsec "eqSection.*" if(issection(app)){sh.color(2)}
    sh.show(0)

}


// data to compute path from each section to the soma
objectvar eqrvp			
objectvar eqslsoma		
objectvar eqslroot

{create eqSection[100]}

factEqRa = 1//0.298
RaMERGINGMETHOD = 0

proc createEqCell(){local nCluster,appLun,sec,col,i,j,jj,secApp
    
    clusterListApp = $o1
    nCluster = clusterListApp.count
    eqsecpathri = new Vector(nCluster,0.0)		//set/reset
    eqsecpathri0 = new Vector(nCluster,0.0)		//set/reset
    orMinpathri = new Vector(nCluster,10^6)		//set/reset
    orMaxpathri = new Vector(nCluster,0.0)		//set/reset
    eqRa = new Vector(nCluster,0.0)		//set/reset
        
    TotOrSurf = orSurf.sum(1,orSurf.size-1)
    TotEqSurf = eqSurf.sum(1,eqSurf.size-1)


    for (i = 0; i < nCluster; i = i + 1){
        
           for jj=0,clusterListApp.o(i).size()-1{
                secApp = clusterListApp.o(i).x(jj)
                eqRa.set(i,eqRa.x(i)+secRa.x(secApp))

                if (secpathri.x(secApp)>orMaxpathri.x(i)){
                   orMaxpathri.set(i,secpathri.x(secApp))
                }
                if (secpathri0.x(secApp)<orMinpathri.x(i)){
                   orMinpathri.set(i,secpathri0.x(secApp))
                }                
           }
           if (clusterListApp.o(i).size()==1){
              if (parsec.x(secApp)!=secApp){
                 orMinpathri.set(i,secpathri.x(parsec.x(secApp)))
              }else{orMinpathri.set(i,0)}
              orMaxpathri.set(i,secpathri.x(secApp))
           }
           
        if(optPrint) print "orMinpathri di ",i,"è: ", orMinpathri.x(i)
        if(optPrint) print "orMaxpathri di ",i,"è: ", orMaxpathri.x(i)
      
 
        eqRa.set(i,eqRa.x(i)/clusterListApp.o(i).size())
        

        
        for (j = 1; j < nCluster; j = j + 1){
            if (clusterParent.x(j)==i) {eqSection[i] connect eqSection[j](0), clusterParentPos.x(j)} 
        }       
        aaa = 0
        if (SURFFACT_PARTIAL_TOTAL==0 && i>0) {surfFact.set(i,TotOrSurf/TotEqSurf) surfFactRmCm.set(i,TotOrSurf/TotEqSurf)}
        eqSection[i] {L = eqL.x(i) diam = eqdiam.x(i) 
                        if($2==0){RaEq = factEqRa*eqRa.x(i)}\
                        else if($2==1){RaEq = factEqRa*eqri2.x(i)*PI*(diam/2)^2*100/L}\
                        else if($2==2){RaEq = factEqRa*(orMaxpathri.x(i)-orMinpathri.x(i))*PI*(diam/2)^2*100/L}
                        /*psection()*/ }
        if (i>=0){      
           eqSection[i] {insert pas e_pas=Vrest \
           g_pas=factRm*surfFactRmCm.x(i)*g_pas0 Ra=RaEq cm=factCm*cm0*surfFactRmCm.x(i)}
        }else if(i==0 || i==1){
           eqSection[i] {insert pas e_pas=Vrest \
           g_pas=surfFactRmCm.x(i)*g_pas0 Ra=RaEq cm=surfFactRmCm.x(i)*cm0}        
        }


    }
    
    totEq=0
    forsec "eqSection.*" {
           nseg = int((L/(d_lambda*lambda_f(freqEq))+0.9)/2)*2 + 1 
           totEq=totEq+nseg
    }   
    if(optPrint) print "totSegEq: ", totEq
  
    sec = 1
    forsec "eqSection.*" if (!issection("eqSection[0]")){
	  {eqslsoma=new SectionList() eqrvp=new RangeVarPlot("v")
	   {eqSection[0]  eqrvp.begin(.5)}  eqrvp.end(.5) eqrvp.list(eqslsoma) }
	  {eqslroot=new SectionList() eqrvp=new RangeVarPlot("v")
	   {ss=0 forsec eqslsoma{ {if(ss==1) eqrvp.begin(.5)}  ss=ss+1}}
	   eqrvp.end(.5)  eqrvp.list(eqslroot)}
	   
	  nsecPath = 0 
	  forsec eqslroot {nsecPath=nsecPath+1 for(x) if(x>0) eqsecpathri.set(sec,eqsecpathri.x(sec)+ri(x))}
      nsecPath1 = 0
      forsec eqslroot {
             nsecPath1=nsecPath1+1 
             for(x) if(x>0 && nsecPath1<nsecPath) eqsecpathri0.set(sec,eqsecpathri0.x(sec)+ri(x))
      }      
        
	  sec = sec+1
    }    
    forsec "eqSection[0]" {for(x) if(x>0) eqsecpathri.set(0,eqsecpathri.x(0)+ri(x))}

    appLun = eqSection[0].L
    {access eqSection[0]}
    {pt3dclear()}
    {pt3dadd(300,0,0,eqSection[0].diam)}
    {pt3dadd(300,appLun,0,eqSection[0].diam)}
    appLun = eqSection[1].L
    {access eqSection[1]}
    {pt3dclear()}
    {pt3dadd(300,0,0,eqSection[1].diam)}
    {pt3dadd(300+appLun,0,0,eqSection[1].diam)}
    appLun = eqSection[2].L
    {access eqSection[2]}
    {pt3dclear()}
    {pt3dadd(300,0,0,eqSection[2].diam)}
    {pt3dadd(300,appLun,0,eqSection[2].diam)}
    
    vbox0 = new VBox()  
    vbox1 = new VBox()  
    strdef strApp 
    
    if(optPrint2) {
        vbox1.intercept(1)
        shPlot = new PlotShape()
        shPlot.scale(-70, 0)  //-70,30,-50
        shPlot.exec_menu("Show Diam")
        shPlot.exec_menu("Shape Plot")
        vbox1.intercept(0)
        vbox1.map(strApp, 400, 295, 500, 420)
    }

    if(optPrint1) {
        vbox0.intercept(1)
        sh = new Shape()
        vbox0.intercept(0)
        sprint(strApp, "Plotting %s", str)
        vbox0.map(strApp, 0, 630, 300, 180)
      
        app = ""
        coloringClusterList(clusterListApp)
        sh.action("app = secname() secApp=-1 forall {secApp=secApp+1 if (issection(app)) break} print secApp,app coloringClusterList(clusterListApp)")
        sh.show(0)  
    }
}


objref ff,ffMorph
objectvar synSecVect, synPosVect, synGmaxVect
objectvar synSecMorphVect, synPosMorphVect, synGmaxMorphVect

proc onfile(){
     nCluster = $o1.count
     ff = new File("./outEq.txt") 
     ff.wopen("./outEq.txt")
     ffMorph = new File("./outMorph.txt") 
     ffMorph.wopen("./outMorph.txt")
     ff.printf("%s\n",strFile)
     ffMorph.printf("%s\n",strFile)
     ff.printf("%d\n",nCluster)
     ind = 0
     forsec "eqSection.*" { 
        ff.printf("%d\t%d\t%-8.8g\t%-8.8g\t%-8.8g\t%-8.8g\t%-8.8g\n",ind,nseg,L,Ra,diam,cm,g_pas)
        ind=ind+1
     }   
     
     ff.printf("%d\t%d\t%d\n",noise0,poissonSeed,positionSeed) 
     ffMorph.printf("%d\t%d\t%d\n",noise0,poissonSeed,positionSeed)      
 
     
     ff.printf("%d\n",synNum) 
     ffMorph.printf("%d\n",synNum) 
     
     //print synNum
     
     if (synNum>0){
         for(i = 0; i < synNum; i = i + 1){
             ff.printf("%d\t%-8.8g\t%-8.8g\n",synSecVect.x(i),synPosVect.x(i),synGmaxVect.x(i))
             ffMorph.printf("%d\t%-8.8g\t%-8.8g\n",synSecMorphVect.x(i),synPosMorphVect.x(i),synGmaxMorphVect.x(i))
         }
     }


     for (i = 0; i < nCluster; i = i + 1){
         ff.printf("%-8.8g\n",surfFact.x(i))
     }     
     
     ffMorph.printf("%-6.6g\n",freq)

     
     ff.close("./outEq.txt")
     ffMorph.close("./outMorph.txt")

}


/*---------------------------------------------
                   ACCURACY FUNCTION
  ---------------------------------------------*/
  
accu = -1
fact2 = 0.35
TmaxNeg = 10
TmaxPos = 10

// Useful variables
objectvar vecTTSpikeOrApp, vecTTSpikeEqApp

func accuracy(){local i, j, nn
     
    tstop = t 
    vecTTSpikeOrApp = $o1
    vecTTSpikeEqApp = $o2
    
    
	TP = 0
	TN = 0
	FP = 0
	FN = 0
	
	if(vecTTSpikeOrApp.size()==0) vecTTSpikeOrApp.append(tstop)
	if(vecTTSpikeEqApp.size()==0) vecTTSpikeEqApp.append(tstop)
	
	for i = 0,vecTTSpikeOrApp.size()-1{
        if (i == 0){tprec = 0
        }else{tprec = vecTTSpikeOrApp.x(i - 1)}
        if (i == vecTTSpikeOrApp.size()-1){tsucc = tstop
        }else{tsucc = vecTTSpikeOrApp.x(i + 1)}
        if (fact2*(vecTTSpikeOrApp.x(i) - tprec) <= TmaxPos){DeltaT1 = fact2*(vecTTSpikeOrApp.x(i) - tprec)
        }else{DeltaT1 = TmaxPos}
        if (i == 0){initNeg = 0
        }else{initNeg = tprec + DeltaT1}
        fineNeg = vecTTSpikeOrApp.x(i) - DeltaT1
        initPos = fineNeg
        if (fact2*(tsucc - vecTTSpikeOrApp.x(i)) <= TmaxPos){DeltaT2 = fact2*(tsucc - vecTTSpikeOrApp.x(i))
        }else{DeltaT2 = TmaxPos}        
        finePos = vecTTSpikeOrApp.x(i) + DeltaT2
        
        for (nn = initNeg; nn <= fineNeg; nn = nn + TmaxNeg) {
            
            inizioInter = nn
			if (nn + TmaxNeg <= fineNeg) {fineInter = nn + TmaxNeg
            }else(fineInter = fineNeg)	
            
            nEq = 0
            for j = 0,vecTTSpikeEqApp.size()-1{
                if (vecTTSpikeEqApp.x(j) > inizioInter && vecTTSpikeEqApp.x(j) < fineInter) nEq = nEq+1
            }
            if (nEq == 0) {TN = TN + 1
            }else{FP = FP + nEq}   
        }	

        nEq = 0
        for j = 0,vecTTSpikeEqApp.size()-1{
            if (vecTTSpikeEqApp.x(j) >= initPos && vecTTSpikeEqApp.x(j) <= finePos) nEq = nEq+1
        }
        
        
        if (nEq == 0) { FN = FN + 1
        }else{TP = TP + 1 FP = FP + nEq - 1}        	


    }
		
    initNeg = finePos
	fineNeg = tstop
     
    
    for (nn = initNeg; nn <= fineNeg; nn = nn + TmaxNeg) {
        
        inizioInter = nn
		if (nn + TmaxNeg <= fineNeg) {fineInter = nn + TmaxNeg
        }else(fineInter = fineNeg)	
        
        nEq = 0
        for j = 0,vecTTSpikeEqApp.size()-1{
            if (vecTTSpikeEqApp.x(j) > inizioInter && vecTTSpikeEqApp.x(j) < fineInter) nEq = nEq+1
        }
        if (nEq == 0) {TN = TN + 1
        }else{FP = FP + nEq}   
    }    
   
    if(TP + TN + FP + FN != 0){ accu = (TP + TN)/(TP + TN + FP + FN) }
        	
    
    return accu
}


surfFactVect = new Vector(100,1.0)		//set/reset

Loading data, please wait...