/* Author: Xu Zhang @UConn, Jul., 2018
Motor cortical neurons: PYN and FSI (abbreviated as PY and FS in the script, respectively)
Modified from Destexhe et al., 1998
*/

// To execute this script individually, please specify the following variables, e.g.:
// PY2PYrnd = 106 // rng seed
// PYrnd = 107
// PYNintrnd = 207
// FSIintrnd = 307
// noiseSwitch=1
// dt = 0.0125

celsius = 36

// Low-threshold spiking pyramidal neurons * 20

create PYcell[20]
objref ocPY[20]

for i = 0,19 {
	access PYcell[i]

	PYcell[i].Ra = 100		// geometry 
	PYcell[i].nseg = 1
	PYcell[i].diam = 96
	PYcell[i].L = 96			// such that area is about 58000 um2
	PYcell[i].cm = 1

	insert pas
	insert mchh2
	insert mcIm
	insert mcCad
	insert mcIt

	// pas
	// PYcell[i].e_pas = -70
	// PYcell[i].g_pas = 1e-5
	PYcell[i].ek = -100
	PYcell[i].ena = 50
	// hh2
	PYcell[i].vtraub_mchh2 = -55	// Resting Vm, BJ was -55
	// PYcell[i].gnabar_hh2 = 0.05	// McCormick=15 muS, thal was 0.09
	// PYcell[i].gkbar_hh2 = 0.005	// spike duration of pyr cells
	// im
	taumax_im = 1000
	// PYcell[i].gkbar_im = 3e-5		// specific to LTS pyr cell
	// cad
	PYcell[i].depth_mcCad = 1		// McCormick= 0.1 um
	PYcell[i].taur_mcCad = 5		// McCormick=1 ms !!!
	PYcell[i].cainf_mcCad = 2.4e-4	// McCormick=0
	PYcell[i].kt_mcCad = 0		// no pump
	// it
	PYcell[i].cai = 2.4e-4 
	PYcell[i].cao = 2 
	PYcell[i].eca = 120 
	// PYcell[i].gcabar_it = 0.001	// specific to LTS pyr cell

	// Change parameters to LTS
	PYcell[i].gcabar_mcIt = 4e-4
	PYcell[i].gkbar_mcIm = 3e-5
	PYcell[i].g_pas = 1e-5
	PYcell[i].e_pas = -85
	PYcell[i].gnabar_mchh2 = 0.05
	PYcell[i].gkbar_mchh2 = 0.005

	PYcell[i] ocPY[i] = new IClamp(0.5)
	ocPY[i].amp = 0.17 // (nA) Simulate baseline activity
	ocPY[i].del = 0
	ocPY[i].dur = 1e10
}

// Connect PY2PY
objref PY2PYsyn[40]

// Random interconnections among PYNs
objref rngPYNint
rngPYNint = new Random(PYNintrnd)
rngPYNint.normal(1,0.04)

// Synaptic noise within PY
objref PYi_noisyn
PYi_noisyn = new Random(PY2PYrnd)
PYi_noisyn.normal(0, 1e-6*noiseSwitch)

PYcell[0] PY2PYsyn[0] = new tanhSyn(0.6)
setpointer PY2PYsyn[0].vpre, PYcell[19].v(0.5)
PY2PYsyn[0].g = PY2PYsyn[0].g * abs(rngPYNint.repick()) * 100

PYcell[0] PY2PYsyn[20] = new tanhSyn(0.6)
setpointer PY2PYsyn[20].vpre, PYcell[1].v(0.5)
PY2PYsyn[20].g = PY2PYsyn[20].g * abs(rngPYNint.repick()) * 100

for i = 1,18 {
	PYcell[i] PY2PYsyn[i] = new tanhSyn(0.6)
	setpointer PY2PYsyn[i].vpre, PYcell[i-1].v(0.5)
	PY2PYsyn[i].g = PY2PYsyn[i].g * abs(rngPYNint.repick()) * 100
	
	PYcell[i] PY2PYsyn[i+20] = new tanhSyn(0.6)
	setpointer PY2PYsyn[i+20].vpre, PYcell[i+1].v(0.5)
	PY2PYsyn[i+20].g = PY2PYsyn[i+20].g * abs(rngPYNint.repick()) * 100
}

PYcell[19] PY2PYsyn[19] = new tanhSyn(0.6)
setpointer PY2PYsyn[19].vpre, PYcell[18].v(0.5)
PY2PYsyn[19].g = PY2PYsyn[19].g * abs(rngPYNint.repick()) * 100

PYcell[19] PY2PYsyn[39] = new tanhSyn(0.6)
setpointer PY2PYsyn[39].vpre, PYcell[0].v(0.5)
PY2PYsyn[39].g = PY2PYsyn[39].g * abs(rngPYNint.repick()) * 100

for i = 0,19 {
	// Synaptic noise
	PYi_noisyn.play(&PY2PYsyn[i].noise)
}

// Fast-spiking interneurons * 2

create FScell[2]
objref OCFS[2]

for i = 0,1 {
	access FScell[i]

	FScell[i].Ra = 100		// geometry 
	FScell[i].nseg = 1
	FScell[i].diam = 67
	FScell[i].L = 67			// such that area is about 29000 um2
	FScell[i].cm = 1

	insert pas
	insert mchh2

	// pas
	FScell[i].ek = -100
	FScell[i].ena = 50
	// hh2
	FScell[i].vtraub_mchh2 = -55	// Resting Vm, BJ was -55

	FScell[i].g_pas = 0.00015
	FScell[i].e_pas = -70
	FScell[i].gnabar_mchh2 = 0.05
	FScell[i].gkbar_mchh2 = 0.01

	FScell[i] OCFS[i] = new IClamp(0.5)
	OCFS[i].amp = 0.15 // (nA) Simulate background excitation from prefrontal cortex
	OCFS[i].del = 0
	OCFS[i].dur = 1e10
}

// Connect FS2FS
objref FS2FSsyn[2]
for i = 0,1 {
	FScell[i] FS2FSsyn[i] = new tanhSyn(0.6)
	FS2FSsyn[i].alpha = 5
	FS2FSsyn[i].tau = 5.56
	FS2FSsyn[i].e = -80
	FS2FSsyn[i].g = 0.0015
}

setpointer FS2FSsyn[0].vpre, FScell[1].v(0.5)
setpointer FS2FSsyn[1].vpre, FScell[0].v(0.5)

for i = 0,1 {
	// Synaptic noise
	PYi_noisyn.play(&FS2FSsyn[i].noise)
}

// Connect PY2FS
objref PY2FSsyn[40]
for i = 0,19 {
	FScell[0] PY2FSsyn[i] = new tanhSyn(0.7+0.01*i)
	setpointer PY2FSsyn[i].vpre, PYcell[i].v(0.5)
	PY2FSsyn[i].g = 0.0025
	
	FScell[1] PY2FSsyn[i+20] = new tanhSyn(0.7+0.01*i)
	setpointer PY2FSsyn[i+20].vpre, PYcell[i].v(0.5)
	PY2FSsyn[i+20].g = 0.0025
}

// Connect FS2PY
objref FS2PYsyn[40]

// Random interconnections FSI-PYN
objref rngFSIint
rngFSIint = new Random(FSIintrnd)
rngFSIint.normal(1,0.0256)

for i = 0,19 {
	PYcell[i] FS2PYsyn[i] = new tanhSyn(0.7)
	setpointer FS2PYsyn[i].vpre, FScell[0].v(0.5)
	FS2PYsyn[i].alpha = 5
	FS2PYsyn[i].tau = 5.56
	FS2PYsyn[i].e = -80
	FS2PYsyn[i].g = 2.2e-2 * abs(rngFSIint.repick()) // * 0.8
	
	PYcell[i] FS2PYsyn[i+20] = new tanhSyn(0.8)
	setpointer FS2PYsyn[i+20].vpre, FScell[1].v(0.5)
	FS2PYsyn[i+20].alpha = 5
	FS2PYsyn[i+20].tau = 5.56
	FS2PYsyn[i+20].e = -80
	FS2PYsyn[i+20].g = 2.2e-2 * abs(rngFSIint.repick()) // * 0.8
}


for i = 0,39 {
	// Synaptic noise
	PYi_noisyn.play(&PY2FSsyn[i].noise)
	PYi_noisyn.play(&FS2PYsyn[i].noise)
}

// Membrane noise in PY
objref MC_noisc, PYnoise[20], FSnoise[2]
MC_noisc = new Random(PYrnd)
MC_noisc.normal(0, 1e-1*noiseSwitch)
for i = 0,19 {
	PYcell[i] PYnoise[i] = new NoisyCurrent(0.2)
	MC_noisc.play(&PYnoise[i].noise)
}
// Membrane noise in FS
for i = 0,1 {
	FScell[i] FSnoise[i] = new NoisyCurrent(0.2)
	MC_noisc.play(&FSnoise[i].noise)
}