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

	CURRENT-CLAMP SIMULATIONS OF TC CELLS
	=====================================
	Accession number ModelDB: 232876

	Simulations of a simplified model of thalamic relay cell.

	The current model is described in the following paper:
		Zeldenrust F, Chameau P, Wadman WJ. Spike and Burst Coding in Thalamocortical Relay Cells. PLoS Computational Biology, 2017.

	The original model, the basis of the current one, is described in the following paper:
	   Destexhe A, Neubig M, Ulrich D and Huguenard JR.  Dendritic
	   low-threshold calcium currents in thalamic relay cells.  
	   Journal of Neuroscience 18: 3574-3588, 1998.

	Please cite these references when using this model.
----------------------------------------------------------------------------*/

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

load_file("nrngui.hoc")		// updated command version of above
nrncontrolmenu()		// create control menu

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 addshape() { local ii	// define subroutine to add a new shape
				// addshape()
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new PlotShape()
	g[ii].scale(-130,50)
}


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

trans = 00

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


//----------------------------------------------------------------------------
//  create multi-compartment geometry and insert currents
//----------------------------------------------------------------------------

xopen("cells/tc3.geo")		// read geometry file

corrD = 7.954			// dendritic surface correction
				// (estimated by fitting voltage-clamp traces)
				
corrR = 0.75  // make R 25% higher, so g=1/1.25=0.8, do also for C for correction RC time
corrC = 0.9   // make RC time smaller

G_pas = corrR*3.79e-5
E_pas = -73			// to fit current-clamp data (was -71 to -73)
E_pas = -76.5			// within 3 mV error
C_m   = corrC*corrR*0.88


forall { 				// insert passive current everywhere
	insert pas
	g_pas = G_pas * corrD
	e_pas = E_pas
	cm = C_m * corrD
	Ra = 173
	L = L
	
}

soma {
	g_pas = G_pas
	cm = C_m

	insert hh2			// insert fast spikes
	ena = 50
	ek = -100
	vtraub_hh2 = -52
	gnabar_hh2 = 0.1
	gkbar_hh2 = 0.1
	
	insert iL     		// insert IL (high-threshold calcium)
	pcabar_iL=0.00005
	
	
	
	insert iC   // IKC (calcium-dependent potassium)
	gkbar_iC=0.001
	}


forall {
	insert itGHK		// T-current everywhere
	cai = 2.4e-4 
	cao = 2 
	eca = 120 
	shift_itGHK = -1	// screening charge shift + 3 mV error
	gcabar_itGHK = corrD * 0.0002
	qm_itGHK = 2.5
	qh_itGHK = 2.5
	

	insert cad			// calcium diffusion everywhere
	depth_cad = 0.1 * corrD
	kt_cad = 0			// no pump
	kd_cad = 1e-4
	taur_cad = 5
	cainf_cad = 2.4e-4	
	
	insert iar 			// Ih everywhere
	ghbar_iar = 1.3e-04

}



xopen("loc3.oc")		// load procedures for localizing T-current

// high distal density, match detailed model in voltage-clamp
localize(1.7e-5,corrD*9.5e-5)

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

Dt = 0.1
npoints = 9030000         // longer traces  

dt = 0.1			// must be submultiple of Dt
tstart = trans
tstop = trans + npoints * Dt
runStopAt = tstop
steps_per_ms = 1/Dt

celsius=30

v_init = -74			// approximate resting Vm

//----------------------------------------------------------------------------
//  insert electrodes in the soma
//----------------------------------------------------------------------------

if(ismenu==0) {
  xopen("El.oc")		// Electrode with series resistance
  ismenu = 1
}

access soma

// DC current electrode for setting holding potential
//----------------------------------------------------------------------------

objectvar El1			// insert electrode
El1 = new Electrode()
electrodes_present = 1


// CURRENT-CLAMP MODE
soma El1.stim.loc(0.5)		// put electrode in current-clamp mode
El1.stim.del = 0
El1.stim.dur = tstop - tstart
El1.stim.amp = 0

//----------------------------------------------------------------------------
//  insert high conductance state
//----------------------------------------------------------------------------

access soma

objref fl
fl = new Gfluct2(0.5)
fl.E_i = -80	

// turn off high-conductance state
//----------------------------------------------------------------------------	
//fl.g_e0    = 0
//fl.g_i0    = 0
//fl.std_e   = 0
//fl.std_i   = 0

// inhibitory high-conductance state
//----------------------------------------------------------------------------	
//fl.g_e0    = 0
//fl.g_i0    = 0.07
//fl.std_e   = 0
//fl.std_i   = 0.01

//create a seed
SEED_0 = 178205489
SEED_1 =127949444
fl.new_seed(SEED_0)

//----------------------------------------------------------------------------
//  insert Poisson spike train(s)
//----------------------------------------------------------------------------

random_stream_offset_ = 1000 
// template RandomStream, modeldB entry 83319
//----------------------------------------------------------------------------
begintemplate RandomStream

public r, repick, start, stream
external random_stream_offset_
objref r
proc init() {
	stream = $1
	r = new Random()
	start()
}
func start() {
	return r.MCellRan4(stream*random_stream_offset_ + 1)
}
func repick() {
	return r.repick()
}
endtemplate RandomStream



MAXSTREAM=1000
random_stream_offset_ = MAXSTREAM

objref Espikesource, rs// our label for excitatory source
dend1[0]  Espikesource = new NetStim(0.5) // arbitrary location

rs = new RandomStream(0)                //random stream
Espikesource.noiseFromRandom(rs.r)      //associate random stream and Espikesource
rs.r.negexp(1) // must specify negexp distribution with mean = 1

Espikesource.interval =100  // average inter-spike interval of the input (in ms), so 10 Hz
Espikesource.number =((tstop - tstart)*1000)/Espikesource.interval
Espikesource.start = 10
Espikesource.noise = 1 // =0 periodic input, =1 Poisson
Espikesource.seed(SEED_0)

// Excitatory synapse on the neuron.
//----------------------------------------------------------------------------
objref Esynapse
dend1[0]  Esynapse = new ExpSyn(0.5) // the location of synapse on the dendrite
Esynapse.tau = 3 // synaptic time constant (ms)
Esynapse.e = 0 // synaptic reversal potential (mV)
 
//connect source and target
objref Econn
thresh = 10 // not important when connection is from NetStim
delay = 35

// strength exponential synapse
// without high conductance state, neuron starts responding to Poisson input from about weight 0.007, at about 0.025 it responds to almost all epscs
// with the high conductance state, it responds to about 80% at 0.1, starts responding at about 0.05

Eweight= 0.1// connection strength in µS
//Eweight = 0

Econn = new NetCon(Espikesource, Esynapse, thresh, delay, Eweight)   

// Inhibitory synapse on the neuron.
//----------------------------------------------------------------------------
objref Esynapse_i
dend1[0]  Esynapse_i = new ExpSyn(0.5) // the location of synapse on the dendrite
Esynapse_i.tau = 10 // synaptic time constant (ms)
Esynapse_i.e = -80 // synaptic reversal potential (mV)
 
//connect source and target
objref Econn_i
thresh_i = 10 // not important when connection is from NetStim
delay_i = 0

// strength inhibitory synapse
Eweight_i= 0.3// connection strength in µS
//Eweight_i = 0

Econn_i = new NetCon(Espikesource, Esynapse_i, thresh_i, delay_i, Eweight_i) 


//prodecure to reset the random class
 proc restart() { 
  rs.start()
}


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

addgraph("soma.v(0.5)",-100,50)			// Vm soma	
addgraph("Esynapse.i",-1e-2,1e-2) 	    // Excitatory Poisson input current	
addgraph("fl.i",-1e-2,1e-2) 			// fluctuating synapse (high conductance state)


//----------------------------------------------------------------------------
//  Save
//----------------------------------------------------------------------------
objref rect, recv, recinternalcalcium, recicap, recih, recit, recikc, recihighcond, recmit, rechit, reco1ih, reco2ih, recmikc, idvec, tvec, recisynapse, recgsynapse, recisynapse_i, recgsynapse_i  

// General
//----------------------------------------------------------------------------
//rect = new Vector()
//rect.record(&t)
recv = new Vector()
recv.record(&soma.v(0.5))
recinternalcalcium = new Vector()
recinternalcalcium.record(&soma.cai(0.5))

// Currents
//----------------------------------------------------------------------------
// recicap = new Vector()
// recicap.record(&soma.i_cap(0.5))
// recih = new Vector()
// recih.record(&soma.ih(0.5))
recit = new Vector()
recit.record(&soma.ica(0.5))
recikc = new Vector()
recikc.record(&soma.ik(0.5))
recihighcond = new Vector()
recihighcond.record(&fl.i)

// Poisson spike train
//----------------------------------------------------------------------------
tvec = new Vector()
idvec = new Vector()
idvec.add(1) 
idvec.mark(g, tvec, "|", 10, 2, 1)
Econn.record(tvec, idvec)
recisynapse = new Vector()
recisynapse.record(&Esynapse.i)
recisynapse_i = new Vector()
recisynapse_i.record(&Esynapse_i.i)


// Gating variables
//----------------------------------------------------------------------------
recmit = new Vector()
recmit.record(&soma.m_itGHK(0.5))
rechit = new Vector()
rechit.record(&soma.h_itGHK(0.5))
// reco1ih = new Vector()
// reco1ih.record(&soma.o1_iar(0.5))
// reco2ih = new Vector()
// reco2ih.record(&soma.o2_iar(0.5))
// recmikc = new Vector()
// recmikc.record(&soma.m_iC(0.5))

// Conductances
//----------------------------------------------------------------------------
recgsynapse = new Vector()
recgsynapse.record(&Esynapse.g)
recgsynapse_i = new Vector()
recgsynapse_i.record(&Esynapse_i.g)


// Procedures for saving
//----------------------------------------------------------------------------
objref savdata
objref tempmatrix, tempmatrix2
tempmatrix = new Matrix()
tempmatrix2 = new Matrix()

proc fileopen() {
    savdata = new File()

    savdata.wopen("tcrelay3c_high_conductance_Esyn_ie.dat")
}
    
proc slaop() {
    tempmatrix.resize(recv.size(),11)   
    tempmatrix.setcol(0,recv)
    tempmatrix.setcol(1,recinternalcalcium)
    tempmatrix.setcol(2,recit)
    tempmatrix.setcol(3,recikc)
    tempmatrix.setcol(4,recihighcond)
    tempmatrix.setcol(5, recisynapse)
    tempmatrix.setcol(6, recgsynapse)
    tempmatrix.setcol(7, recisynapse_i)
    tempmatrix.setcol(8, recgsynapse_i)
    tempmatrix.setcol(9, recmit)
    tempmatrix.setcol(10, rechit)
	tempmatrix.fprint(savdata, " %g")

    savdata.close()
}