// 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 April 2015
//

begintemplate TTC

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

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

objref synlist, SecSyn, ASCIIrpt, HDF5rpt, APC
objref all, somatic, apical, basal, dends
objref preconlist,synLocList, rList,preTrainList
objref rd1
objref this
objref sref,fih
objref locVec

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()
  dends = 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()

  apic[0] dends.append()
  apic[1] dends.append()
  dend[0] dends.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

  soma {nseg = 1 Ra = 100 diam = 13.47 L = 23.17 cm = 1 V = -80 e_pas = -90}
  dend {nseg = 1 Ra = 100 diam = 10.28 L = 282.13 cm = 2 V = -80 e_pas = -90}
  apic[0] {nseg = 1 Ra = 100 diam = 5.04 L = 700.0 cm = 2 V = -80 e_pas = -90}
  apic[1] {nseg = 1 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()
  synLocList = new List()
  stimThreshold =0
  preTrainList = 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
}

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

double siteVec[2]

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

	lengthA = apic[0].L + apic[1].L
	lengthB = dend.L
        //lengthA = 0
        //lengthB = 0
        //forsec "apical" {
        //  lengthA = lengthA + L
        //}
        //forsec "basal" {
        //  lengthB = lengthB + L
        //}
        pA = lengthA/(lengthA + lengthB)


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

		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]
		locVec = new Vector()
		locVec.append(treeInd)
		locVec.append(siteVec[0])
		locVec.append(siteVec[1])
		synLocList.append(locVec)

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

		sprint(cmd2,"%s[siteVec[0]] sref = new SectionRef()",treename)
		execute(cmd2,this)

		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
				synlist.o[syni].weight_factor_NMDA = weight_factor_NMDA
	
				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
			}
		} 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 = 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
	strdef cmd2,treename
	
	TconVec = $o1

	for(ii=0;ii<TconVec.size();ii+=1){
		if(TconVec.x[ii]!=0){
			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 vectors
proc setpretrains(){local j
  for(j=0;j<(NsynE+NsynI);j+=1){
    preTrainList.append($o1.o[j])
  }
}

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

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

endtemplate TTC