/*
sj3-cortex.hoc

Defined in this file:
Synapses
Cell Templates
Exogenous Feeds
----------------------------------------
Cells are instantiated
Basic connector defined and dipole set-up
*/

  create acell_home_
  access acell_home_
L=1
diam=1

/////////////////////////////////////
// Receptor Mechanisms realized with Exp2Syn. These can be replaced 
// with kinetic schemes
/////////////////////////////////////
begintemplate AMPA
public s
objref s
proc init() {
s = new Exp2Syn(0.5)
s.tau1=0.5
s.tau2=5
s.e=0
}
endtemplate AMPA

begintemplate NMDA
public s
objref s
proc init() {
s = new Exp2Syn(0.5)
s.tau1=1
s.tau2=20
s.e=0
}
endtemplate NMDA

begintemplate GABAA
public s
objref s
proc init() {
s = new Exp2Syn(0.5)
s.tau1=0.5
s.tau2=5
s.e=-80
}
endtemplate GABAA

begintemplate GABAB
public s
objref s
proc init() {
s = new Exp2Syn(0.5)
s.tau1=1
s.tau2=20
s.e=-80
}
endtemplate GABAB

///////////////////////////////////////
// Template for Layer V Pyranidal Cell
///////////////////////////////////////
begintemplate Layer5_pyr
public is_art
public init, topol, basic_shape, subsets, geom, biophys
public x, y, z, position, connect2target
public ampa,nmda,gabaa,gabab

public soma, dend
public all, somatic, dendritic, apical, basal, dend0, dend1, dend2, dend3
public dend4, dend5, dend6, dend7
// if gbar_x is declared as 'public', can not be set in biophys() below
// public gbar_ar, gbar_cat 

objref synlist

proc init() {
  topol()
  subsets()
  geom()
  biophys()
  geom_nseg()
  synlist = new List()
  synapses()
  x = y = z = 0 // only change via position
}

create soma, dend[8]

proc topol() { local i
  connect dend(0), soma(1)
  for i = 1, 2 connect dend[i](0), dend(1)
  for i = 3, 4 connect dend[i](0), dend[i-1](1)
  connect dend[5](0), soma(0) //was soma(1)this is correct! 
  for i = 6, 7 connect dend[i](0), dend[5](1)
  basic_shape()
}

 proc basic_shape() {
// THESE AND LENGHTHS MUST CHANGE TOGETHER!!!
  soma {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(0, 23, 0, 1)}
  dend {pt3dclear() pt3dadd(0, 23, 0, 1) pt3dadd(0, 83, 0, 1)}
  dend[1] {pt3dclear() pt3dadd(0, 83, 0, 1) pt3dadd(-150, 83, 0, 1)}
  dend[2] {pt3dclear() pt3dadd(0, 83, 0, 1) pt3dadd(0, 483, 0, 1)}
  dend[3] {pt3dclear() pt3dadd(0, 483, 0, 1) pt3dadd(0, 883, 0, 1)}
  dend[4] {pt3dclear() pt3dadd(0, 883, 0, 1) pt3dadd(0, 1133, 0, 1)}
  dend[5] {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(0, -50, 0, 1)}
  dend[6] {pt3dclear() pt3dadd(0, -50, 0, 1) pt3dadd(-106, -156, 0, 1)}
  dend[7] {pt3dclear() pt3dadd(0, -50, 0, 1) pt3dadd(106, -156, 0, 1)}
} 

objref all, somatic, dendritic, apical, basal, dend0, dend1, dend2, dend3, dend4
objref dend5, dend6, dend7
proc subsets() { local i
  objref all, somatic, dendritic, apical, basal, dend0, dend1, dend2, dend3 
  objref dend4, dend5, dend6, dend7

  all = new SectionList()
    soma all.append()
    for i=0, 7 dend[i] all.append()

  somatic = new SectionList()
    soma somatic.append()

  dendritic = new SectionList()
    for i=0, 7 dend[i] dendritic.append()

  apical = new SectionList()
    for i=0, 4 dend[i] apical.append()

  basal = new SectionList()
    for i=5, 7 dend[i] basal.append()

  dend0 = new SectionList()
	dend[0] dend0.append()

  dend1 = new SectionList()
	dend[1] dend1.append()

  dend2 = new SectionList()
	dend[2] dend2.append()

  dend3 = new SectionList()
	dend[3] dend3.append()

  dend4 = new SectionList()
	dend[4] dend4.append()

  dend5 = new SectionList()
	dend[5] dend5.append()

  dend6 = new SectionList()
	dend[6] dend6.append()

  dend7 = new SectionList()
	dend[7] dend7.append()


}
//increased lengths and diams by 70% for human size
proc geom() {
  forsec all {  }
    // soma.L = 13 //BUSH 1999 spike amp smaller
   soma.L=39 //Bush 1993
   dend.L = 102
   dend[1].L = 255
   dend[2].L = 680 //default 400
   dend[3].L = 680 //default 400
   dend[4].L = 425
   dend[5].L = 85
   dend[6].L = 255 // default 150
   dend[7].L = 255 // default 150
    //soma.diam = 18.95 //Bush 1999
   soma.diam = 28.9 //Bush 1993
   dend.diam = 10.2
   dend[1].diam = 5.1
   dend[2].diam = 7.48 //default 4.4
   dend[3].diam = 4.93 //default 2.9
   dend[4].diam = 3.4
   dend[5].diam = 6.8
   dend[6].diam = 8.5
   dend[7].diam = 8.5
}
proc geom_nseg() {
  soma area(.5) // make sure diam reflects 3d points
   forsec all { if (L < 50) {nseg=1} else {nseg=int(L/50)} }
}
proc biophys() {
/* USING DEFAULT NEURON HH GIVES CORRECT SPIKE WIDTH */
	forsec all{
		Ra=200
		cm=0.85 // decreased by 70%
		
	}
	forsec somatic {
		insert hh
		gnabar_hh=0.16 // to match latency with old segmentation
		gkbar_hh=0.01 //0.01 0.09 for spiking
		gl_hh=0.0000426 //decreased by 70%
		el_hh=-65 //default -65
		insert ca
		gbar_ca=60// 100 //150 pS/um2=0.00015 S/cm2
		insert cad
		taur_cad=20 //20 ms 
		insert kca
		gbar_kca=0.0002//0.0002 //(0.00015 mho/cm2= S/cm2)
		insert km
		gbar_km=200//200 //10 pS/um2
                insert cat
                gbar_cat=0.002
                insert ar
                gbar_ar=0.00006
	}
	forsec dendritic {
		insert hh
		gnabar_hh=0.14 //0.45 S/cm2
		gkbar_hh=0.01 //0.09 for spiking
		gl_hh=0.0000426
		el_hh=-71 //default -71
		insert ca
		gbar_ca=60//100 //15 pS/um2
		insert cad
		taur_cad=20 //20 ms 
		insert kca
		gbar_kca=0.0002//0.0002 //0.015(mho/cm2= S/cm2)
		insert km
		gbar_km=200//200 //10 pS/um2
                insert cat
                gbar_cat=0
                insert ar
                gbar_ar=0

	} 

 }  //end proc biophysics 
proc position() { local i
  soma for i = 0, n3d()-1 {
    pt3dchange(i, $1-x+x3d(i), $2-y+y3d(i), $3-z+z3d(i), diam3d(i))
  }
  x = $1  y = $2  z = $3
}
proc connect2target() { //$o1 target point process, $o2 returned NetCon
  soma $o2 = new NetCon(&v(1), $o1)
}

objref ampa[9],nmda[9],gabaa[9],gabab[9]
proc synapses() {
  dend[0]{ ampa[0] = new AMPA() nmda[0] = new NMDA() gabaa[0] = new GABAA() gabab[0] = new GABAB() }
  dend[1]{ ampa[1] = new AMPA() nmda[1] = new NMDA() gabaa[1] = new GABAA() gabab[1] = new GABAB() }           
  dend[2]{ ampa[2] = new AMPA() nmda[2] = new NMDA() gabaa[2] = new GABAA() gabab[2] = new GABAB() }           
  dend[3]{ ampa[3] = new AMPA() nmda[3] = new NMDA() gabaa[3] = new GABAA() gabab[3] = new GABAB() }           
  dend[4]{ ampa[4] = new AMPA() nmda[4] = new NMDA() gabaa[4] = new GABAA() gabab[4] = new GABAB() }           
  dend[5]{ ampa[5] = new AMPA() nmda[5] = new NMDA() gabaa[5] = new GABAA() gabab[5] = new GABAB() }           
  dend[6]{ ampa[6] = new AMPA() nmda[6] = new NMDA() gabaa[6] = new GABAA() gabab[6] = new GABAB() }           
  dend[7]{ ampa[7] = new AMPA() nmda[7] = new NMDA() gabaa[7] = new GABAA() gabab[7] = new GABAB() }           
  soma   { ampa[8] = new AMPA() nmda[8] = new NMDA() gabaa[8] = new GABAA() gabab[8] = new GABAB() }
}

func is_art() { return 0 }

endtemplate Layer5_pyr

//////////////////////////////////////////////
/////////////////////////////////////////////
//Template for Layer II Pyramidal Cell
//////////////////////////////////////////////

begintemplate Layer2_pyr
public init, topol, basic_shape, subsets, geom, biophys
public synlist, x, y, z, position, connect2target
public ampa,nmda,gabaa,gabab

public soma, dend
public all, somatic, dendritic, apical, basal, dend0, dend1
public dend2, dend3, dend4, dend5, dend6 
public pre, ampa, gabaa, gabab, nmda
public gbar_ar, gbar_cat

objref synlist

proc init() {
  topol()
  subsets()
  geom()
  biophys()
  geom_nseg()
  synlist = new List()
  synapses()
  x = y = z = 0 // only change via position
}

create soma, dend[7]

proc topol() { local i
  connect dend(0), soma(1)
  for i = 1, 2 connect dend[i](0), dend(1)
  connect dend[3](0), dend[2](1)
  connect dend[4](0), soma(0) //was soma(1), 0 is correct!
  for i = 5, 6 connect dend[i](0), dend[4](1)
  basic_shape()
}
proc basic_shape() {
  soma {pt3dclear() pt3dadd(-50, 765, 0, 1) pt3dadd(-50, 778, 0, 1)}
  dend {pt3dclear() pt3dadd(-50, 778, 0, 1) pt3dadd(-50, 813, 0, 1)}
  dend[1] {pt3dclear() pt3dadd(-50, 813, 0, 1) pt3dadd(-250, 813, 0, 1)}
  dend[2] {pt3dclear() pt3dadd(-50, 813, 0, 1) pt3dadd(-50, 993, 0, 1)}
  dend[3] {pt3dclear() pt3dadd(-50, 993, 0, 1) pt3dadd(-50, 1133, 0, 1)}
  dend[4] {pt3dclear() pt3dadd(-50, 765, 0, 1) pt3dadd(-50, 715, 0, 1)}
  dend[5] {pt3dclear() pt3dadd(-50, 715, 0, 1) pt3dadd(-156, 609, 0, 1)}
  dend[6] {pt3dclear() pt3dadd(-50, 715, 0, 1) pt3dadd(56, 609, 0, 1)}
}

objref all, somatic, dendritic, apical, basal, dend0, dend1
objref dend2, dend3, dend4, dend5, dend6
proc subsets() { local i
  objref all, somatic, dendritic, apical, basal, dend0, dend1
  objref dend2, dend3, dend4, dend5, dend6
  all = new SectionList()
    soma all.append()
    for i=0, 6 dend[i] all.append()

  somatic = new SectionList()
    soma somatic.append()

  dendritic = new SectionList()
    for i=0, 6 dend[i] dendritic.append()

  apical = new SectionList()
    for i=0, 3 dend[i] apical.append()

  basal = new SectionList()
    for i=4, 6 dend[i] basal.append()

  dend0 = new SectionList()
    dend dend0.append()

  dend1 = new SectionList()
    dend[1] dend1.append()

  dend2 = new SectionList()
    dend[2] dend2.append()

  dend3 = new SectionList()
    dend[3] dend3.append()

  dend4 = new SectionList()
    dend[4] dend4.append()

  dend5 = new SectionList()
    dend[5] dend5.append()

  dend6 = new SectionList()
    dend[6] dend6.append()

}
// increased by 70% for human
proc geom() {
  forsec all {  }
   soma.L = 22.1
   dend.L = 59.5
   dend[1].L = 340
   dend[2].L = 306
   dend[3].L = 238
   dend[4].L = 85
   dend[5].L = 255
   dend[6].L = 255
   
    soma.diam = 23.4
    dend.diam = 4.25
    dend[1].diam = 3.91
    dend[2].diam = 4.08
    dend[3].diam = 3.4
    dend[4].diam = 4.25
    dend[5].diam = 2.72
    dend[6].diam = 2.72

}
proc geom_nseg() {
  soma area(.5) // make sure diam reflects 3d points
   forsec all { if (L < 50) {nseg=1} else {nseg=int(L/50)} }
}
proc biophys() {
/* USING DEFAULT NEURON HH GIVES CORRECT SPIKE WIDTH */
 forsec all {
    Ra = 200
    cm = 0.6195 // decreased 70% from 2.065
	}
forsec somatic {
		insert hh
		gnabar_hh=0.18 // to match latency with old segmentation
		gkbar_hh=0.01 //0.09
		gl_hh=0.0000426 //decreased 70% for human
		el_hh=-65
		insert km
		gbar_km=250 //450 pS/um2
                insert cat
                gbar_cat=0.0
                insert ar
                gbar_ar=0.0

	}
	forsec dendritic {
		insert hh
		gnabar_hh=0.15//0.45 S/cm2
		gkbar_hh=0.01 //0.09
		gl_hh=0.0000426 
		el_hh=-65
		insert km
		gbar_km=250 //450 pS/um2 =uS/cm2
                insert cat
                gbar_cat=0.0
                insert ar
                gbar_ar=0.0

	} 

 }  //end proc biophysics 


proc position() { local i
  soma for i = 0, n3d()-1 {
    pt3dchange(i, $1-x+x3d(i), $2-y+y3d(i), $3-z+z3d(i), diam3d(i))
  }
  x = $1  y = $2  z = $3
}
proc connect2target() { //$o1 target point process, $o2 returned NetCon
  soma $o2 = new NetCon(&v(1), $o1)
}

objref ampa[8],nmda[8],gabaa[8],gabab[8]
proc synapses() {
  dend[0]{ ampa[0] = new AMPA() nmda[0] = new NMDA() gabaa[0] = new GABAA() gabab[0] = new GABAB() }
  dend[1]{ ampa[1] = new AMPA() nmda[1] = new NMDA() gabaa[1] = new GABAA() gabab[1] = new GABAB() }
  dend[2]{ ampa[2] = new AMPA() nmda[2] = new NMDA() gabaa[2] = new GABAA() gabab[2] = new GABAB() }
  dend[3]{ ampa[3] = new AMPA() nmda[3] = new NMDA() gabaa[3] = new GABAA() gabab[3] = new GABAB() }
  dend[4]{ ampa[4] = new AMPA() nmda[4] = new NMDA() gabaa[4] = new GABAA() gabab[4] = new GABAB() }
  dend[5]{ ampa[5] = new AMPA() nmda[5] = new NMDA() gabaa[5] = new GABAA() gabab[5] = new GABAB() }
  dend[6]{ ampa[6] = new AMPA() nmda[6] = new NMDA() gabaa[6] = new GABAA() gabab[6] = new GABAB() }
  soma   { ampa[7] = new AMPA() nmda[7] = new NMDA() gabaa[7] = new GABAA() gabab[7] = new GABAB() }
} 

func is_art() { return 0 }

endtemplate Layer2_pyr

/////////////////////////////////////
// Template for Exogenous feeds to network
begintemplate FeedX
public pp, connect2target, x, y, z, position, is_art
external acell_home_
objref pp

proc init() {
  acell_home_ pp = new NetStimG(.5)
    pp.MeanInterval = 10
    pp.SD = 0.5
    pp.MeanStart = 25
    pp.StartSD = 1
    pp.number = 10
}
func is_art() { return 1 }
proc connect2target() { $o2 = new NetCon(pp, $o1) }
proc position(){x=$1  y=$2  z=$3}
endtemplate FeedX

/////////////////////////////////////
// Template for Inhibitory cells
/////////////////////////////////////
begintemplate Inhib
public cell,connect2target,init,ampa,nmda,gabaa,gabab
external AMPA,NMDA,GABAA,GABAB

objref ampa[1],nmda[1],gabaa[1],gabab[1]
create cell

proc init(){
  access  cell
                L=39
                diam=20
                Ra=200
                cm=0.85
                insert hh
synapses()
}

proc synapses() {
cell {
      ampa = new AMPA()
      nmda = new NMDA()
      gabaa = new GABAA()
      gabab = new GABAB()
     }
}

// proc connect2target() {
// $o2 = new NetCon(cell(0.5), $o1)
proc connect2target() { //$o1 target point process, $o2 returned NetCon
  cell $o2 = new NetCon(&v(1), $o1)
}

endtemplate Inhib
/////////////////////////////////////
/////////////////////////////////////
// Keep track of connections
objref    AMPAconnects, NMDAconnects, GABAAconnects, GABABconnects

AMPAconnects = new List()
NMDAconnects = new List()
GABAAconnects = new List()
GABABconnects = new List()
objref netcon

/////////////////////////////////////
///////The Network Constructor/////////////
// Called by procedures defined in wiring_proc.hoc
//////////////////////////////////////////
proc rconnect(){
// usage:       rconnect(source-cell, target-cell, target-section, rtype,
// weight, delay, threshold )
// rtype= 0=AMPA ||| 1=NMDA ||| 2=GABAa ||| 3= GABAb
// -----------------------------------------------------------------------

if ( $4 == 0){ $o1.connect2target($o2.ampa[$3].s, netcon)
netcon.weight = $5 netcon.delay = $6 netcon.threshold = $7
 AMPAconnects.append(netcon) }

if ( $4 == 1){ $o1.connect2target($o2.nmda[$3].s, netcon)
netcon.weight = $5 netcon.delay = $6 netcon.threshold = $7
 NMDAconnects.append(netcon) }

if ( $4 == 2){ $o1.connect2target($o2.gabaa[$3].s, netcon)
netcon.weight = $5 netcon.delay = $6 netcon.threshold = $7
 GABAAconnects.append(netcon) }

if ( $4 == 3){ $o1.connect2target($o2.gabab[$3].s, netcon)
netcon.weight = $5 netcon.delay = $6 netcon.threshold = $7
 GABABconnects.append(netcon) }

if (( $4 > 3) || ($3 < 0)){ print "Receptor Type not supported" }
}

////////////////////////////////////////
// Create some cells
////////////////////////////////////////

objref PL5[10], IPL5[3], PL2[10], IPL2[3]
// Create 10 Layer V Pyramidal Cells
for i=0,9 { PL5[i] = new Layer5_pyr() PL5[i].position(i*100,0,0) } 
//  and 3 Layer V inhibitory
for i=0,2 { IPL5[i] = new Inhib() }     
// Create 10 Layer II/III  Pyramidal Cells
for i=0,9 { PL2[i] = new Layer2_pyr() PL2[i].position(i*100,0,i*100) } 
//  and 3 Layer II/III inhibitory
for i=0,2 { IPL2[i] = new Inhib() }     


objref  FF, FB, FF2
FF = new FeedX() // Thalamic input (Feed-forward) 
 FF.pp.MeanInterval = 1
 FF.pp.SD=0
 FF.pp.MeanStart = 25
 FF.pp.StartSD = 2.5
 FF.pp.number = 1

FB = new FeedX() // Feed-Back (eg. Pre-frontal input) 
 FB.pp.MeanInterval = 1
 FB.pp.SD=0
 FB.pp.MeanStart = 70
 FB.pp.StartSD = 6
 FB.pp.number = 1

FF2 = new FeedX() // Thalamic input (Feed-forward) 
 FF2.pp.MeanInterval = 1
 FF2.pp.SD=0
 FF2.pp.MeanStart = 135
 FF2.pp.StartSD = 7
 FF2.pp.number=1

///////////////////////////////////////////
// Procedures for keeping lists of connections categorized by receptor type
///////////////////////////////////////////
objref strobj
strobj = new StringFunctions()
objref vbox
proc gabaa_browser(){
GABAAconnects.browser("GABAa")
GABAAconnects.select_action("Param_Window_GABAA(hoc_ac_)")
}

proc Param_Window_GABAA(){
vbox = new VBox()
vbox.intercept(1)
print "presynaptic:"
strobj.references(GABAAconnects.o($1).pre)
print "postsynaptic:"
strobj.references(GABAAconnects.o($1).postcell)

xpanel("GABAa $1")
xvarlabel("strobj.references(GABAAconnects.o($1).postcell)")
xpvalue("weight",&GABAAconnects.o($1).weight,1,"")
xpvalue("delay",&GABAAconnects.o($1).delay,1,"")
xpvalue("threshold",&GABAAconnects.o($1).threshold,1,"")

xpanel()
vbox.intercept(0)
vbox.map
}

proc gabab_browser(){
GABABconnects.browser()
GABABconnects.select_action("Param_Window_GABAB(hoc_ac_)")
}

proc Param_Window_GABAB(){
vbox = new VBox()
vbox.intercept(1)
print "presynaptic:"
strobj.references(GABABconnects.o($1).pre)
print "postsynaptic:"
strobj.references(GABAAconnects.o($1).postcell)

xpanel("GABAb $1")
xvarlabel("")
xpvalue("weight",&GABABconnects.o($1).weight,1,"")
xpvalue("delay",&GABABconnects.o($1).delay,1,"")
xpvalue("threshold",&GABABconnects.o($1).threshold,1,"")

xpanel()
vbox.intercept(0)
vbox.map
}

proc ampa_browser(){
AMPAconnects.browser()
AMPAconnects.select_action("Param_Window_AMPA(hoc_ac_)")
}

proc Param_Window_AMPA(){
vbox = new VBox()
vbox.intercept(1)
print "presynaptic:"
strobj.references(AMPAconnects.o($1).pre)
print "postsynaptic:"
strobj.references(AMPAconnects.o($1).postcell)

xpanel("AMPA $1")
xvarlabel("")
xpvalue("weight",&AMPAconnects.o($1).weight,1,"")
xpvalue("delay",&AMPAconnects.o($1).delay,1,"")
xpvalue("threshold",&AMPAconnects.o($1).threshold,1,"")

xpanel()
vbox.intercept(0)
vbox.map
}
proc nmda_browser(){
NMDAconnects.browser()
NMDAconnects.select_action("Param_Window_NMDA(hoc_ac_)")
}

proc Param_Window_NMDA(){
vbox = new VBox()
vbox.intercept(1)
print "presynaptic:"
strobj.references(NMDAconnects.o($1).pre)
print "postsynaptic:"
strobj.references(NMDAconnects.o($1).postcell)

xpanel("NMDA $1")
xvarlabel("")
xpvalue("weight",&NMDAconnects.o($1).weight,1,"")
xpvalue("delay",&NMDAconnects.o($1).delay,1,"")
xpvalue("threshold",&NMDAconnects.o($1).threshold,1,"")

xpanel()
vbox.intercept(0)
vbox.map
}

/////////////////////////////////////
// DIPOLE procedures
/////////////////////////////////////
load_file("dipole.hoc")


objref dipoleL5, dipoleL2
dipoleL5 = new List()
dipoleL2 = new List()

for i=0, 9 {PL5[i].soma dipoleL5.append(new DipoleForCell() )  }
for i=0, 9 {PL2[i].soma dipoleL2.append(new DipoleForCell() )  }

dL5=0
func L5_dipole(){ local dp
dp=0
for i=0, 9{ forsec PL5[i].all {dp += Qsum_dipole} }
dL5 = dp
return dp
}
 
dL2=0
func L2_dipole(){ local dp
dp=0
for i=0, 9{ forsec PL2[i].all {dp += Qsum_dipole} }
dL2=dp
return dp
}

func dset() {local s
        s=0
        forsec $o1 {s+=Qsum_dipole}
        return s
}
/////////////////////////////////////////////////////////////////////////