/*******************************************************************
Load a cell's morphology file and then load Typc C biophysical model
2/18/2005
*******************************************************************/
/*--------------------playing random seed----------------*/
use_mcell_ran4(1)
lowindex = 1001
mcell_ran4_init(lowindex)
highindex = 1

print "------------------------------Type C Model-----------------------"

//flag for wether or not to use NMDA-R: 0 for without NMDA-R, 1 for with NMDA-R
//this flag will be used in subsequent hoc files: NregNsyn.hoc, Nregsyn.hoc etc
NMDAFLAG = 0

if (NMDAFLAG == 0) { //only AMPA-R
    wtAmpa = 1
    wtNmda = 0
} else { //NMDAFLAG = 1, with both AMPA-R, NMDA-R
    wtAmpa = 0.8
    wtNmda = 0.2
}

xopen("./c20466.hoc")             // geometry file

soma area(0.5) //make sure diam reflects 3d points


//USER5MAX, APIDENDMAX, BASDENDMAX are set in cell's hoc file
TOTALDEND = USER5MAX + APIDENDMAX + BASDENDMAX + 3 //"3": counting start at 0
objref dend[TOTALDEND]


/****************************************************************************
Grouping dendritic sections together in an array dend, later for being used to
randomly put synapses on the proximal apical, distal apical and basal dendrite
*****************************************************************************/
index = 0

forsec "dendrite" { //this will includes both dendrites and apical-dendrites
    dend[index] = new SectionRef()
    index += 1
}

forsec "user5" {
    dend[index] = new SectionRef()
    index += 1
} //now index = USER5MAX + 1

//quadratic function fitting result A and B are set in cell's hoc file
A = A1
B = B1

//synapse rise and decay time constants
TAU1=0.2 //ms
TAU2=10  //ms

//cvode_active(1)

dist=1
rel=0.1

Rm = 28000  //unit: ohm-cm^2
RmDend = Rm/1
RmSoma = Rm
RmAx = Rm

Cm    = 1
CmSoma= Cm
CmAx  = Cm
CmDend = Cm*1

RaAll= 50
RaSoma=50  
RaAx = 50

Vrest = -65
gna =  .02
AXONM = 2
gkdr = 0.01
celsius = 35.0  
KMULT =  0.025
KMULTP = 0.025
gcan=0.0//005
gcal=0.0//005
gcat=0.0//005
ghd=0.00005
nash=0


forsec "axon" {insert pas e_pas=Vrest g_pas = 1/RmAx Ra=RaAx cm=CmAx}
forsec "soma" {insert pas e_pas=Vrest g_pas = 1/RmSoma Ra=RaSoma cm=CmSoma}
forsec "dendrite"{insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}
forsec "user5" {insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}

access soma
freq=50
load_file("fixnseg.hoc")           
geom_nseg()

tot=0
forall {tot=tot+nseg}
distance()
maxdist=0
forsec "user5" for(x) {if (distance(x)>maxdist) {maxdist=distance(x)}}
print "total # of segments (50Hz): ",tot, "  max path distance: ", maxdist

objref  sref, blist[USER5MAX+1],   aplist
strdef dend2, trunk

for i=0, USER5MAX {
    blist[i] = new Vector()
}
aplist = new Vector(APIDENDMAX+1)


//*********mapping bifurcations******************


for i=0, APIDENDMAX apical_dendrite[i] {
	while (!issection("user5.*")) {
//	    print "before  ", i, secname()
	    sref = new SectionRef()
	    access sref.parent
	    sprint(dend2, secname())
	}
//	print "apical ",i," ",dend2
	for k=0, USER5MAX user5[k] {
	    sprint(trunk,secname())
	    x=strcmp(dend2, trunk)
	    if (x==0) {blist[k].append(i) aplist.x[i]=k}
	}
}
//************************************************


forsec "axon" {   
                insert nax gbar_nax=gna * AXONM	sh_nax=nash
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP*0.2
}

forsec "soma" {   
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 ar_na3=1 sh_na3=nash gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
}

forsec "dendrite" {
	    if(!issection("apical_dendrite.*")) {
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 ar_na3=1 sh_na3=nash gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP

	    }
}

forsec "user5" {

        insert ds
        insert hd ghdbar_hd=ghd
        insert na3 ar_na3=1 gbar_na3=gna sh_na3=nash
        insert kdr gkdrbar_kdr=gkdr
        insert kap gkabar_kap=0
        insert kad gkabar_kad=0

        for (x) if (x>0 && x<1) {
                xdist = distance(x)
                if (xdist>350) {xdist=350}
                ghdbar_hd(x) = ghd*(1+3*xdist/100)
                        if (xdist > 100){
                                vhalfl_hd=-81
                                gkabar_kad(x) = KMULT*(1+xdist/100)
                                } else {
                                vhalfl_hd=-73
                                gkabar_kap(x) = KMULTP*(1+xdist/100)
                        }
//		print secname()," ",x, gkabar_kap(x), gkabar_kad(x)
        }
}


for i=0, APIDENDMAX apical_dendrite[i] {
        insert ds
        insert hd
        insert na3 ar_na3=1 gbar_na3=gna sh_na3=nash
        insert kdr gkdrbar_kdr=gkdr
        insert kap
        insert kad

        gkabar_kad = user5[aplist.x[i]].gkabar_kad(1)
        gkabar_kap = user5[aplist.x[i]].gkabar_kap(1)
        vhalfl_hd = user5[aplist.x[i]].vhalfl_hd
        ghdbar_hd = user5[aplist.x[i]].ghdbar_hd(1)

//      print secname()," ", gkabar_kap, gkabar_kad, ghdbar_hd
}

//correct the default section name
access soma

proc init() {
	t=0
        forall {
        v=Vrest
        if (ismembrane("nax") || ismembrane("na3")) {ena=55}
        if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {ek=-90}
        if (ismembrane("hd") ) {ehd_hd=-30}
	}
	finitialize(Vrest)
        fcurrent()

        forall {
	for (x) {
	if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)}
	if (ismembrane("hd")) {e_pas(x)=e_pas(x)+i_hd(x)/g_pas(x)}
		}
	}
	cvode.re_init()
	cvode.event(tstop)
	access soma
//	g.begin()
}


proc advance() {
	fadvance()
//	g.plot(t)
//	g.flush()
//	p.flush()
//	doNotify()
}

init()