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

	Simplified model of regular-spiking cortical neuron
	===================================================

        Single-compartment model of "regular-spiking" pyramidal neurons,
        which is the most commonly encountered electrophysiological type
        of excitatory cell in cortex.  The model is based on the presence
	of three voltage-dependent currents: 
        - INa, IK: action potentials
        - IM: slow K+ current for spike-frequency adaptation
        (no ICa/IK[Ca] in this model)

  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.

  The model was taken from a thalamocortical model, described in:

   Destexhe, A., Contreras, D. and Steriade, M.
   Mechanisms underlying the synchronizing action of corticothalamic
   feedback through inhibition of thalamic relay cells.
   J. Neurophysiol. 79: 999-1016, 1998.



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

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


batch = 0		// for batch processing (use nrniv)


//----------------------------------------------------------------------------
//  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("sPY_template")		// read geometry file

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

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









//----------------------------------------------------------------------------
//  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.75		// ** was 0.9

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 = 300
	El[i].stim.dur = 400
	El[i].stim.amp = CURR_AMP
}


electrodes_present=1





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

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

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






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

strdef gtxt

if(batch == 0) {
  addgraph("PY[0].soma.m_im",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)
}