/*----------------------------------------------------------------------------

	Simplified model of bursting cortical neuron
	============================================

        Single-compartment model of bursting pyramidal neurons.
	The model is based on the presence of four
        voltage-dependent currents: 
        - INa, IK: action potentials
        - IM: slow K+ current for spike-frequency adaptation
        - IL: L-type calcium currents for burst generation

	Repetitive bursting with stronger IL

  Model described in:

   Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., 
   Bal, T., Fregnac, Y., Markram, H. and Destexhe, A.
   Minimal Hodgkin-Huxley type models for different classes of
   cortical and thalamic neurons.
   Biological Cybernetics 99: 427-441, 2008.


        Alain Destexhe, CNRS, 2009
	http://cns.iaf.cnrs-gif.fr

----------------------------------------------------------------------------*/


//----------------------------------------------------------------------------
//  load and define general graphical procedures
//----------------------------------------------------------------------------

load_file("stdrun.hoc")

objectvar g[20]			// max 20 graphs
ngraph = 0

proc addgraph() { local ii	// define subroutine to add a new graph
				// addgraph("variable", minvalue, maxvalue)
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph()
	g[ii].size(tstart,tstop,$2,$3)
	g[ii].xaxis()
	g[ii].yaxis()
	g[ii].addvar($s1,1,0)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
}

proc addtext() { local ii	// define subroutine to add a text graph
				// addtext("text")
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph()
	g[ii].size(0,tstop,0,1)
	g[ii].xaxis(3)
	g[ii].yaxis(3)
	g[ii].label(0.1,0.8,$s1)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
	text_id = ii
}

proc addline() {		// to add a comment to the text window
				// addline("text")
	g[text_id].label($s1)
}


if(ismenu==0) {
  nrnmainmenu()			// create main menu
  nrncontrolmenu()		// crate control menu
}


//----------------------------------------------------------------------------
//  transient time
//----------------------------------------------------------------------------

trans = 0000

print " "
print ">> Transient time of ",trans," ms"
print " "









//----------------------------------------------------------------------------
//  create PY cells
//----------------------------------------------------------------------------

print " "
print "<<==================================>>"
print "<<            CREATE CELLS          >>"
print "<<==================================>>"
print " "

xopen("sPYbr_template")		// read geometry file

ncells = 1			// nb of cells in each layer <<>>

objectvar PY[ncells]
for i=0,ncells-1 {
  PY[i] = new sPYbr()
}









//----------------------------------------------------------------------------
//  Rebound parameters
//----------------------------------------------------------------------------

// L-current density:
// 0.0001 -> RS behavior
// 0.00017 -> one burst, then RS
// 0.0002 -> two bursts, then RS
// 0.00022 -> repetitive bursting

PY[0].soma {
	gcabar_ical = 2.2e-4
	gkbar_im = 3e-5
	g_pas = 1e-5
	e_pas = -85
	gnabar_hh2 = 0.05
	gkbar_hh2 = 0.005
}






//----------------------------------------------------------------------------
//  insert electrode in each PY cell
//----------------------------------------------------------------------------

if(ismenu==0) {
  load_file("electrod.hoc")	// electrode template
  ismenu = 1
}

objectvar El[ncells]			// create electrodes

CURR_AMP = 0.15

for i=0,ncells-1 {			// insert one in each cell
	PY[i].soma El[i] = new Electrode()
	PY[i].soma El[i].stim.loc(0.5)
	El[i].stim.del = 500
	El[i].stim.dur = 2000
	El[i].stim.amp = CURR_AMP
}


electrodes_present=1




//----------------------------------------------------------------------------
//  setup simulation parameters
//----------------------------------------------------------------------------

Dt = .1				// macroscopic time step <<>>
npoints = 30000

dt = 0.1			// must be submultiple of Dt
tstart = trans
tstop = trans + npoints * Dt
runStopAt = tstop
steps_per_ms = 5
celsius = 36
v_init = -84






//----------------------------------------------------------------------------
//  add graphs
//----------------------------------------------------------------------------

strdef gtxt

if(batch == 0) {
  addgraph("PY[0].soma.m_im",0,1)
//  addgraph("PY[0].soma.m_iahp",0,1)
  for i=0,ncells-1 {
	sprint(gtxt,"PY[%d].soma.v(0.5)",i)
	addgraph(gtxt,-120,40)
  }
}





//----------------------------------------------------------------------------
//  add text
//----------------------------------------------------------------------------

access PY[0].soma

proc text() {
  sprint(gtxt,"%d PY cells",ncells)
  addtext(gtxt)
  sprint(gtxt,"Passive: gleak=%g Eleak=%g",PY.soma.g_pas,PY.soma.e_pas)
  addline(gtxt)
  sprint(gtxt,"HH: gNa=%g, gK=%g, vtraub=%g",PY.soma.gnabar_hh2,\
  PY.soma.gkbar_hh2,PY.soma.vtraub_hh2)
  addline(gtxt)
  sprint(gtxt,"IM: g=%g, taumax=%g",PY.soma.gkbar_im,taumax_im)
  addline(gtxt)
  sprint(gtxt,"Ca++: tau=%g, depth=%g, cainf=%g",taur_cad,depth_cad,cainf_cad)
  addline(gtxt)
  sprint(gtxt,"IL: g=%g",PY.soma.gcabar_ical)
  addline(gtxt)
}