// Author: Etay Hay 2014
// Dendritic excitability and gain control in recurrent cortical microcircuits (Hay and Segev, 2014, Cerebral Cortex)
//
// Cell template of L5 thick-tufted pyramidal cell (TTC)
//
// Modified by Tuomo Maki-Marttunen 2015-2016
//

begintemplate TTC

public init, biophys, geom_nseg, biophys_inhomo
public synlist, connect2target, APC
public locateSites, getLongestBranch, distributeSyn, distributeSyn2
public initRand,insertMCcons
public preconlist ,synLocList, rList, preTrainList, preTrainIndexList
public setnetworkparameters,initPreSynTrain,queuePreTrains,setpretrains
public pA

public soma, dend, apic, stimThreshold, ASCIIrpt, HDF5rpt, getAbsSecIndex
public all, somatic, apical, basal, nSecSoma, nSecApical, nSecBasal, nSecAll, SecSyn
public nsegs
public NsynsE, NsynsI

objref synlist, SecSyn, ASCIIrpt, HDF5rpt, APC
objref all, somatic, apical, basal
objref preconlist, syngroupIndexList, rList, preTrainList, preTrainIndexList
objref rd1, rd2
objref this
objref sref,fih
objref NsynsE, NsynsI
objref synvarVecList, treenameList

strdef tstr

create soma[1], dend[1], apic[2]


proc init() {localobj nl,import
  all = new SectionList()
  somatic = new SectionList()
  basal = new SectionList()
  apical = new SectionList()

  apic[0] apical.append()
  apic[1] apical.append()
  soma[0] somatic.append()
  dend[0] basal.append()

  apic[0] all.append()
  apic[1] all.append()
  soma[0] all.append()
  dend[0] all.append()
  
  connect dend[0](0), soma[0](0)
  connect apic[0](0), soma[0](1)
  connect apic[1](0), apic[0](1)

  access soma

  soma insert pas
  dend insert pas
  apic[0] insert pas
  apic[1] insert pas

  nsegs = 5

  soma {nseg = nsegs Ra = 100 diam = 13.47 L = 23.17 cm = 1 V = -80 e_pas = -90}
  dend {nseg = nsegs Ra = 100 diam = 10.28 L = 282.13 cm = 2 V = -80 e_pas = -90}
  apic[0] {nseg = nsegs Ra = 100 diam = 5.04 L = 700.0 cm = 2 V = -80 e_pas = -90}
  apic[1] {nseg = nsegs Ra = 100 diam = 5.04 L = 600.53 cm = 2 V = -80 e_pas = -90}

  soma {
    insert Ca_LVAst
    insert Ca_HVA
    insert SKv3_1
    insert SK_E2
    insert K_Tst
    insert K_Pst
    insert Nap_Et2
    insert NaTa_t
    insert CaDynamics_E2
    insert Ih
    ek = -85
    ena = 50
    gIhbar_Ih = 0.0002
    g_pas = 0.0000338
    decay_CaDynamics_E2 = 460.0
    gamma_CaDynamics_E2 = 0.000501
    gCa_LVAstbar_Ca_LVAst = 0.00343
    gCa_HVAbar_Ca_HVA = 0.000992
    gSKv3_1bar_SKv3_1 = 0.693
    gSK_E2bar_SK_E2 = 0.0441
    gK_Tstbar_K_Tst = 0.0812
    gK_Pstbar_K_Pst = 0.00223
    gNap_Et2bar_Nap_Et2 = 0.00172
    gNaTa_tbar_NaTa_t = 2.04
  }
  dend {
    insert Ih
    gIhbar_Ih = 0.0002
    g_pas = 0.0000467
  }
  apic[0] {
    insert Ca_LVAst
    insert Ca_HVA
    insert SKv3_1
    insert SK_E2
    insert NaTa_t
    insert Im
    insert CaDynamics_E2
    insert Ih
    ek = -85
    ena = 50
    decay_CaDynamics_E2 = 122
    gamma_CaDynamics_E2 = 0.000509
    gSK_E2bar_SK_E2 = 0.0012
    gSKv3_1bar_SKv3_1 = 0.000261
    gNaTa_tbar_NaTa_t = 0.0213
    gImbar_Im = 0.0000675
    g_pas = 0.0000589
    gIhbar_Ih = 0.0004
    gCa_LVAstbar_Ca_LVAst = 0.0187
    gCa_HVAbar_Ca_HVA = 0.000555
  }
  apic[1] {
    insert Ca_LVAst
    insert Ca_HVA
    insert SKv3_1
    insert SK_E2
    insert NaTa_t
    insert Im
    insert CaDynamics_E2
    insert Ih
    ek = -85
    ena = 50
    decay_CaDynamics_E2 = 122
    gamma_CaDynamics_E2 = 0.000509
    gSK_E2bar_SK_E2 = 0.0012
    gSKv3_1bar_SKv3_1 = 0.000261
    gNaTa_tbar_NaTa_t = 0.0213
    gImbar_Im = 0.0000675
    g_pas = 0.0000589
    gIhbar_Ih = 0.0004
    gCa_LVAstbar_Ca_LVAst = 0.0187
    gCa_HVAbar_Ca_HVA = 0.000555
  }
  geom_nseg()
  area(0.5)
  distance()
  access soma
 
  synlist = new List()
  preconlist = new List()
  stimThreshold =0
  preTrainList = new List()
  preTrainIndexList = new List()
  syngroupIndexList = new List()
  synvarVecList = new List()
  pA = 0.5
}

proc setnetworkparameters(){
  rcpWeightFactor = $1 //how stronger reciprocal connections are on average
  EsynConductance = $2
  IsynConductance = $3
  NsynE = $4
  NsynI = $5
  contactsNumE = $6
  weight_AMPA = $7
  weight_factor_NMDA = $8
  weight_GABA = $9
  gNoiseCoeff = $10
}

proc geom() {
}

proc geom_nseg() {local nSec, L1, L2, D1, D2, nSeg1, nSeg2
  soma area(.5) // make sure diam reflects 3d points
  nSec = 0
  forsec all {
    nseg = 1 + 2*int(L/40)
    nSec = nSec + 1
  }

  nSecAll = 5
  nSecSoma = 1
  nSecApical = 2
  nSecBasal = 1
}

proc connect2target() { //$o1 target point process, $o2 returned NetCon
  soma $o2 = new NetCon(&v(1), $o1)
  $o2.threshold = -30
}

// $s1 section
func getLongestBranch(){local maxL,d localobj distallist,sref
    sprint(tstr,"%s distance()",$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
	}

// $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)    
    
	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
}


//========================================================================================
//================================= microcircuit related stuff============================
//========================================================================================

// $1 some number that is function of the TTC index
//
// Have each TTC with unique sites, but certain portion of inputs identical to root TTC
proc initRand() {
  rList = new List() //for stochastic synapses

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

double siteVec[2]

proc distributeSyn() {local sitenum,i, segInd, compInd localobj sl
	strdef treename,cmd2

        NsynsE = new List()
  	NsynsI = new List()
  	for i = 0, 2 {     
    	  NsynsE.append(new Vector(nsegs))
    	  NsynsI.append(new Vector(nsegs))
  	}

	lengthA = apic[0].L + apic[1].L
	lengthB = dend.L
        pA = lengthA/(lengthA + lengthB)

        for(i=0;i<2*3*nsegs;i+=1) {
          syngroupIndexList.append(new Vector())
        }

	for(i=0;i<(NsynE+NsynI);i+=1){
		if (rd1.repick()<pA){
			treename = "apic"
                        compInd = 1
		} else {
			treename = "dend"
                        compInd = 0
		}

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

  		sitenum = int((sl.count()-1)*rd1.repick())
                compInd = compInd + sl.o[sitenum].x[0] // if we are at apical, and sl.o[sitenum].x[0]=1, then compInd = 2, otherwise 1 at apical, and 0 at basal
                segInd = int(sl.o[sitenum].x[1]*nsegs)             
                if (i<NsynE) {                                                
                  NsynsE.o[compInd].x[segInd] = NsynsE.o[compInd].x[segInd] + 1
  		  syngroupIndexList.o[compInd*nsegs+segInd].append(i)
                } else {                                                       
                  NsynsI.o[compInd].x[segInd] = NsynsI.o[compInd].x[segInd] + 1
  		  syngroupIndexList.o[3*nsegs + compInd*nsegs+segInd].append(i)
                }                                                              
	}
}        


proc distributeSyn2() {local sitenum,syni,preconi,i,i1, segInd, compInd localobj sl,nilstim
	strdef cmd2
        //objref treenameList
        treenameList = new List()

        treenameList.append(new String("dend[0]"))
        treenameList.append(new String("apic[0]"))
        treenameList.append(new String("apic[1]"))

        for(i=0;i<NsynsE.count();i+=1){
            sprint(cmd2,"access %s",treenameList.o[i].s)
            execute(cmd2,this)
            sprint(cmd2,"%s sref = new SectionRef()",treenameList.o[i].s)
            execute(cmd2,this)

            for(iseg=0;iseg<NsynsE.o[i].size();iseg+=1) {
                segx = (0.5+iseg)/NsynsE.o[i].size()
                sref {
                        synlist.append(new ProbAMPANMDA2groupdet(segx))
                        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*gNoiseCoeff
                        synlist.o[syni].weight_factor_NMDA = weight_factor_NMDA

                        Nsyns_thissyn = NsynsE.o[i].x[iseg]
                        synvarVecList.append(new Vector(3*Nsyns_thissyn))
                        for i1=0,Nsyns_thissyn-1 {
                          synvarVecList.o[synvarVecList.count()-1].x(i1) = 0
                          synvarVecList.o[synvarVecList.count()-1].x(Nsyns_thissyn+i1) = 1
                          synvarVecList.o[synvarVecList.count()-1].x(2*Nsyns_thissyn+i1) = 0
                          }
                        synlist.o[syni].setVec(synvarVecList.o[synvarVecList.count()-1])
                        preconlist.append(new NetCon(nilstim, synlist.o[syni]))
                        preconi = preconlist.count()-1 //connection index
                        preconlist.o[preconi].weight = weight_AMPA
                        preconlist.o[preconi].delay = 0
                }
            }
        }
        for(i=0;i<NsynsI.count();i+=1){
            sprint(cmd2,"access %s",treenameList.o[i].s)
            execute(cmd2,this)
            sprint(cmd2,"%s sref = new SectionRef()",treenameList.o[i].s)
            execute(cmd2,this)

            for(iseg=0;iseg<NsynsI.o[i].size();iseg+=1) {
                segx = (0.5+iseg)/NsynsI.o[i].size()
                sref {
                        synlist.append(new ProbUDFsyn2groupdet(segx))
                        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*gNoiseCoeff

                        Nsyns_thissyn = NsynsI.o[i].x[iseg]
                        synvarVecList.append(new Vector(3*Nsyns_thissyn))
                        for i1=0,Nsyns_thissyn-1 {
                          synvarVecList.o[synvarVecList.count()-1].x(i1) = 0
                          synvarVecList.o[synvarVecList.count()-1].x(Nsyns_thissyn+i1) = 1
                          synvarVecList.o[synvarVecList.count()-1].x(2*Nsyns_thissyn+i1) = 0
                          }
                        synlist.o[syni].setVec(synvarVecList.o[synvarVecList.count()-1])
                        preconlist.append(new NetCon(nilstim, synlist.o[syni]))
                        preconi = preconlist.count()-1 //connection index
                        preconlist.o[preconi].weight = weight_GABA
                        preconlist.o[preconi].delay = 0
                }
            }
        }
}


// adds the microcircuit connections
// $o1 conVec - where 0 (no connection), 1 (one way), rcpWeightFactor (reciprocated)
proc insertMCcons(){local sitenum,syni,ii,jj localobj TconVec,sl,NcontVec
	strdef cmd2,treename
	
	TconVec = $o1
        NcontVec = $o2

	for(ii=0;ii<TconVec.size();ii+=1){
		if(TconVec.x[ii]!=0){
                        contactsNumE = NcontVec.x[ii]
			for(jj=0;jj<contactsNumE;jj+=1){
				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,this)
		
				sprint(cmd2,"%s[siteVec[0]] sref = new SectionRef()",treename)
				execute(cmd2,this)
	
				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 = 0
					synlist.o[syni].Fac = 0
					synlist.o[syni].Use = 0.25
					synlist.o[syni].u0 = 0
					synlist.o[syni].gmax = TconVec.x[ii] * EsynConductance
				}
			}
		}
	}
}

//$o1 list of event time vectors
proc setpretrains(){local i,j,k,indsum localobj preTrainThis,preTrainIndicesThis,sortedIndicesThis
  for (i=0;i<syngroupIndexList.count();i+=1) { //The two loops can be combined as syngroupIndexList[0-14] only contain exc. and syngroupIndexList[15-29] only contain inh. syngroup indices
    preTrainThis = new Vector()
    preTrainIndicesThis = new Vector()
    for (j=0;j<syngroupIndexList.o[i].size();j+=1) {
      k = syngroupIndexList.o[i].x[j]
      preTrainThis.append($o1.o[k])
      preTrainIndicesThis.append(new Vector($o1.o[k].size(),j))
    }
    sortedIndicesThis = new Vector()
    sortedIndicesThis = preTrainThis.sortindex()
    preTrainList.append(preTrainThis.index(sortedIndicesThis))    
    preTrainIndexList.append(preTrainIndicesThis.index(sortedIndicesThis))    
    sortedIndicesThis.resize(0) 
    syngroupIndexList.o[i].resize(0)
  }
  syngroupIndexList.remove_all()
}

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

// sets presynaptic spike events
proc initPreSynTrain(){local syni,ti
	for(syni=0;syni<preTrainList.count();syni+=1){
		for(ti=0;ti<preTrainList.o[syni].size();ti+=1){
			preconlist.o[syni].event(preTrainList.o[syni].x[ti])
		}
                synlist.o[syni].setVec2(preTrainIndexList.o[syni])
                //synlist.o[syni].printVec2()
                preTrainList.o[syni].resize(0)
	}
        preTrainList.remove_all()
}

endtemplate TTC