// The axon is the venerable Mainen 1996 axon.  Not perfect, yet very stable.  
// Modified by Tuomo Maki-Marttunen, 2016, basing on TTC.hoc by Etay Hay, Dendritic excitability and gain control in recurrent cortical microcircuits (Hay and Segev, 2014, Cerebral Cortex)

create a_soma

n_axon_seg = 2
create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]

objref preconlist,synLocList, rList, preTrainList
objref rd1
objref sref, fih
objref synlist

proc create_axon() {
  n_axon_seg = 2
  create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]

  a_soma {
    equiv_diam = sqrt(area(.5)/(4*PI))
  }

  iseg {               // initial segment between hillock + myelin
     L = 15
     nseg = 5
     diam = equiv_diam/10        // see Sloper and Powell 1982, Fig.71
  }

  hill {                
    L = 20
    nseg = 5
    diam(0:1) = 2*iseg.diam:iseg.diam
  }

  // construct myelinated axon with nodes of ranvier

  for i=0,n_axon_seg-1 {
    myelin[i] {         // myelin element
      nseg = 2
      L = 100
      diam = iseg.diam
    }
    node[i] {           // nodes of Ranvier
      nseg = 1
      L = 1.0
      diam = iseg.diam*.75       // nodes are thinner than axon
    }
  }

  a_soma connect hill(0), 0.5
  hill connect iseg(0), 1
  iseg connect node[0](0), 1
  node[0] connect myelin[0](0), 1

  for i=0,n_axon_seg-2  { 
      myelin[i] connect node[i+1](0), 1
      node[i+1] connect myelin[i+1](0), 1 
  }
}

// --------
// Spines
// --------
proc add_spines() { local a
  forsec "dend" {
    cm*=$1
    g_pas*=$1
  }
  forsec "apic" {
    cm*=$1
    g_pas*=$1
  }
	
}


proc load_3dcell() {

// $s1 filename

  forall delete_section()
  xopen($s1)
  access a_soma
  
  // make sure no compartments exceed 50 uM length
  forall {
    diam_save = diam
    n = int(L/20)    
    if (n==0) n=1
    nseg = n 
    if (n3d() == 0) diam = diam_save
  }    
  create_axon()
  init_cell()
}


proc init_cell() {

  // passive
  forall {
    insert pas
    Ra = ra 
    cm = c_m 
    g_pas = 1/rm
    e_pas = epas_sim
  }
    // exceptions along the axon
     forsec "myelin" cm = cm_myelin
     forsec "node" g_pas = g_pas_node
  forall {
  	insert iA 
	insert kslow 
	insert na
	insert iH
	insert cah
	insert car
	insert cad
	insert bk
	insert sk
  }

  forall if(ismembrane("k_ion")) ek = Ek
  forall if(ismembrane("na_ion")) ena = Ena    
  forall eh=-33
  forall if(ismembrane("ca_ion")) {
 	vshift_ca = 0
  }
  synlist = new List()
  preconlist = new List()
  preTrainList = new List()

  lengthA = 0
  lengthB = 0
  forsec "apic" {
          lengthA = lengthA + L
  }
  forsec "dend" {
          lengthB = lengthB + L
  }
  pA = lengthA/(lengthA + lengthB)

}

proc density() {

  forall {
    Ra = ra 
    cm = c_m 
    g_pas = 1/rm
    e_pas = epas_sim
    vshiftm_na= na_shift1
    vshifth_na= na_shift2
    taum_scale_na=na_taum_scale
    tauh_scale_na=na_tauh_scale
    q10_iH=ih_q10
    shift_cah=cah_shift
    shifth_cah=cah_shifth

    shift_car=car_shift
    shifth_car=car_shifth
    qm_car=car_qm
  }


  add_spines(2)
   // exceptions along the axon
   forsec "myelin" cm = cm_myelin
   forsec "node" g_pas = g_pas_node
   forall vshift_na = -5
   forsec "myelin" {
  	gbar_na    = gna_soma
	gbar_kslow =gkslow_beta 
	gbar_iA    = gka_beta
   }

  hill {
  	gbar_na = gna_node
	vshiftm_na=7
	vshifth_na=3
	gbar_iA = gka_node
	gbar_kslow = gkslow_node
  }
	
  iseg {
  	gbar_na = gna_node
	vshiftm_na=7
	vshifth_na=3
	gbar_iA = gka_node
	gbar_kslow = gkslow_node
  }
	
  forsec "node" {
  	gbar_na = gna_node
	vshiftm_na=7
	vshifth_na=3
	gbar_iA = gka_node
	gbar_kslow = gkslow_node
  }

   // we put in all of teh dendrites the same values as in the a_soma and then overwrite the appical values	
  forsec "dend" {
    gbar_na = gna_soma
    gbar_kslow = gkslow_start+gkslow_beta
    gbar_iA = gka_start+gka_beta
    gbar_iH=gih_start 
    pbar_car=pcar_soma
    pbar_cah=pcah_soma
    gbar_sk=gsk_dend
    gbar_bk=gbk_dend
  }
  
  a_soma {
    gbar_na = gna_soma
    gbar_kslow = gkslow_start+gkslow_beta
    gbar_iA = gka_start+gka_beta
    gbar_iH=gih_start
    pbar_car=pcar_soma
    pbar_cah=pcah_soma
    gbar_sk=gsk_soma
    gbar_bk=gbk_soma
  }

  
  a_soma distance(0,0.5)
 
  forsec  ApicalDendSectionName {
	dist1=distance(0)
 	dist2=distance(1)

	
	gbar_na = gna_api
	if (dist1<dist_na){
           gbar_na(0:1) = (gna_soma+dist1*(gna_api-gna_soma)/dist_na):(gna_soma+dist2*(gna_api-gna_soma)/dist_na)
        }


	pbar_cah = pcah_api
        if (dist1<dist_cah){
           pbar_cah(0:1) = (pcah_soma+dist1*(pcah_api-pcah_soma)/dist_cah):(pcah_soma+dist2*(pcah_api-pcah_soma)/dist_cah) 
        }

        pbar_car = pcar_api
        if (dist1<dist_car){
           pbar_car(0:1) = (pcar_soma+dist1*(pcar_api-pcar_soma)/dist_car):(pcar_soma+dist2*(pcar_api-pcar_soma)/dist_car) 
	   //(pcar_soma+dist1*(pcar_api-pcar_soma)/dist_car):(pcar_soma+dist2*(pcar_api-pcar_soma)/dist_car)
        } 


	gbar_iA(0:1) = (gka_start+gka_beta*exp(gka_alpha*dist1)):(gka_start+gka_beta*exp(gka_alpha*dist2))

	gbar_kslow(0:1) = (gkslow_start+gkslow_beta*exp(gkslow_alpha*dist1)):(gkslow_start+gkslow_beta*exp(gkslow_alpha*dist2))
	
	gbar_iH(0:1) = (gih_start+gih_end/(1+exp(gih_alpha*(dist1-gih_x2)))):(gih_start+gih_end/(1+exp(gih_alpha*(dist2-gih_x2))))
	 
	gbar_bk = gbk_dend
        if (dist1<dist_bk){
           gbar_bk(0:1) = (gbk_soma+dist1*(gbk_dend-gbk_soma)/dist_bk):(gbk_soma+dist2*(gbk_dend-gbk_soma)/dist_bk)
        }

	gbar_sk = gsk_dend
	if (dist1<dist_sk){
           gbar_sk(0:1) = (gsk_soma+dist1*(gsk_dend-gsk_soma)/dist_sk):(gsk_soma+dist2*(gsk_dend-gsk_soma)/dist_sk)
        }

	
    }

        forsec  ApicalDendSectionName {

                if (pcah_soma<pcah_api){
                        if (pbar_cah>pcah_api){
                                pbar_cah=pcah_api
                        }
                }
                if (pcah_soma>pcah_api){
                        if (pbar_cah<pcah_api){
                                pbar_cah=pcah_api
                        }
                }

                if (pcar_soma<pcar_api){
                        if (pbar_car>pcar_api){
                                pbar_car=pcar_api
                        }
                }
                if (pcar_soma>pcar_api){
                        if (pbar_car<pcar_api){
                                pbar_car=pcar_api
                        }
                }

                if (gbar_na<gna_api) gbar_na=gna_api
		if (gbar_bk<gbk_dend) gbar_bk=gbk_dend
		if (gbar_sk<gsk_dend) gbar_sk=gsk_dend	
        }


}

// Copied from Hay model
// $s1 section
func getLongestBranch(){local maxL,d localobj distallist,sref
  //Commented out: Always calculate the distance with respect to soma 0.5
  //sprint(tstr,"%s distance()",$s1)
  //execute(tstr,this)

  //if(0==strcmp($s1,"axon")){
  //  sprint(tstr,"%s[0] distance(1)",$s1)
  //  execute(tstr,this)
  //}

  maxL = 0
  d = 0
  distallist = new SectionList()
  forsec $s1 {
    sref = new SectionRef()
    if (sref.nchild==0) distallist.append()
  }
  forsec distallist{
    d = distance(1)
    if(maxL<d) maxL = d
  }
  // for the soma case
  if (maxL == 0) {
    $s1 {
      maxL = L
    }
  }
  return maxL
}

//Copied from Hay model
// $s1 section
// $2 distance x in micrometers
// return list of [1,2] vectors  - of the appropriate section and the location in each vector
obfunc locateSites() {local maxL,site,d0,d1,siteX,i localobj vv,ll
  ll = new List()

  //sprint(tstr,"%s distance()",$s1)
  //execute(tstr,this)

  //if(0==strcmp($s1,"axon")){
  //  sprint(tstr,"%s[0] distance(1)",$s1)
  //  execute(tstr,this)
  //}

  maxL = getLongestBranch($s1)
  site = $2
  i = 0
  forsec $s1 {
  if (distance(0) < distance(1)) {
      d0 = distance(0)
      d1 = distance(1)
    } else {
      d1 = distance(0)
      d0 = distance(1)
    }

    if (site <= d1 && site >= d0) {
      siteX = (site-d0)/(d1-d0)
      secNum = i
      vv = new Vector()
      ll.append(vv.append(secNum,siteX))
    }
    i = i+1
  }
  return ll
}


proc setparameters(){
  EsynConductance = $1
  IsynConductance = $2
  NsynE = $3
  NsynI = $4
}

proc initRand() {
  rList = new List() //for stochastic synapses

        rd1 = new Random($1) // unique to this TTC
        rd1.uniform(0,1)
}

double siteVec[2]

proc distributeSyn() {local sitenum,syni,preconi,i localobj sl,nilstim
        strdef treename,cmd2

        for(i=0;i<(NsynE+NsynI);i+=1){
                if (i==(NsynE+NsynI)/4-1) { print "25% of synapses done" }
                if (i==(NsynE+NsynI)/2-1) { print "50% of synapses done" }
                if (i==3*(NsynE+NsynI)/4-1) { print "75% of synapses done" }
                if (rd1.repick()<pA){
                        treename = "apic"
                } else {
                        treename = "dend"
                }

                sl = locateSites(treename,rd1.repick()*getLongestBranch(treename))

                sitenum = int((sl.count()-1)*rd1.repick())
                siteVec[0] = sl.o[sitenum].x[0]
                siteVec[1] = sl.o[sitenum].x[1]

                sprint(cmd2,"access %s[siteVec[0]]",treename)
                execute(cmd2)

		//print "distributesSyn(): ",cmd2
		//print "siteVec = ",siteVec[0],siteVec[1]
                sprint(cmd2,"%s[siteVec[0]] sref = new SectionRef()",treename)
		//print "distributesSyn(): ",cmd2
                execute(cmd2)
                if (i<NsynE){
                        sref {
                                synlist.append(new ProbAMPANMDA2(siteVec[1]))
                                syni = synlist.count()-1 //synapse index
                                rList.append(new Random(int(1000000*rd1.repick())))
                                rList.o[syni].negexp(1)
                                synlist.o[syni].setRNG(rList.o[syni])
                                synlist.o[syni].tau_r_AMPA = 0.3
                                synlist.o[syni].tau_d_AMPA = 3
                                synlist.o[syni].tau_r_NMDA = 2
                                synlist.o[syni].tau_d_NMDA = 65
                                synlist.o[syni].e = 0
                                synlist.o[syni].Dep = 800
                                synlist.o[syni].Fac = 0
                                synlist.o[syni].Use = 0.6
                                synlist.o[syni].u0 = 0
                                synlist.o[syni].gmax = EsynConductance
                                preconlist.append(new NetCon(nilstim, synlist.o[syni]))
                                preconi = preconlist.count()-1 //connection index
                                preconlist.o[preconi].weight = 1
                                preconlist.o[preconi].delay = 0
                        }
                } else {
                        sref {
                                synlist.append(new ProbUDFsyn2(siteVec[1]))
                                syni = synlist.count()-1 //synapse index
                                rList.append(new Random(int(1000000*rd1.repick())))
                                rList.o[syni].negexp(1)
                                synlist.o[syni].setRNG(rList.o[syni])
                                synlist.o[syni].tau_r = 1
                                synlist.o[syni].tau_d = 20
                                synlist.o[syni].e = -80
                                synlist.o[syni].Dep = 800
                                synlist.o[syni].Fac = 0
                                synlist.o[syni].Use = 0.25
                                synlist.o[syni].u0 = 0
                                synlist.o[syni].gmax = IsynConductance
                                preconlist.append(new NetCon(nilstim, synlist.o[syni]))
                                preconi = preconlist.count()-1 //connection index
                                preconlist.o[preconi].weight = 1
                                preconlist.o[preconi].delay = 0
                        }
                }
                if (i==(NsynE+NsynI)-1) { print "synapses done" }
        }
}

proc setsynapseconductances(){local syni, i
        for(i=0;i<(NsynE+NsynI);i+=1){
                if (i<NsynE){
                         synlist.o[i].gmax = $1
                } else {
                         synlist.o[i].gmax = $2
		}
        }
}

//$o1 list of vectors
proc setpretrains(){local j
  for(j=0;j<(NsynE+NsynI);j+=1){
    preTrainList.append($o1.o[j])
  }
}

proc queuePreTrains(){
        fih = new FInitializeHandler("initPreSynTrain()")
}

// sets presynaptic spike events
proc initPreSynTrain(){local ti,si
        for(ti=0;ti<preTrainList.count();ti+=1){
                for(si=0;si<preTrainList.o[ti].size();si+=1){
                        preconlist.o[ti].event(preTrainList.o[ti].x[si])
                }
        }
}