Knox implementation of Destexhe 1998 spike and wave oscillation model (Knox et al 2018)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:234233
" ...The aim of this study was to use an established thalamocortical computer model to determine how T-type calcium channels work in concert with cortical excitability to contribute to pathogenesis and treatment response in CAE. METHODS: The model is comprised of cortical pyramidal, cortical inhibitory, thalamocortical relay, and thalamic reticular single-compartment neurons, implemented with Hodgkin-Huxley model ion channels and connected by AMPA, GABAA , and GABAB synapses. Network behavior was simulated for different combinations of T-type calcium channel conductance, inactivation time, steady state activation/inactivation shift, and cortical GABAA conductance. RESULTS: Decreasing cortical GABAA conductance and increasing T-type calcium channel conductance converted spindle to spike and wave oscillations; smaller changes were required if both were changed in concert. In contrast, left shift of steady state voltage activation/inactivation did not lead to spike and wave oscillations, whereas right shift reduced network propensity for oscillations of any type...."
Reference:
1 . Knox AT, Glauser T, Tenney J, Lytton WW, Holland K (2018) Modeling pathogenesis and treatment response in childhood absence epilepsy. Epilepsia 59:135-145 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex; Thalamus;
Cell Type(s): Thalamus reticular nucleus GABA cell; Thalamus geniculate nucleus/lateral principal GLU cell; Hodgkin-Huxley neuron; Neocortex layer 4 pyramidal cell; Neocortex fast spiking (FS) interneuron;
Channel(s): I h; I Na,t; I K,leak; I T low threshold; I M;
Gap Junctions:
Receptor(s): GabaA; GabaB; AMPA;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Spindles; Oscillations;
Implementer(s): Knox, Andrew [knox at neurology.wisc.edu]; Destexhe, Alain [Destexhe at iaf.cnrs-gif.fr];
Search NeuronDB for information about:  Thalamus geniculate nucleus/lateral principal GLU cell; Thalamus reticular nucleus GABA cell; GabaA; GabaB; AMPA; I Na,t; I T low threshold; I K,leak; I M; I h;
/
KnoxEtAl2017
README.html
README_.txt
ampa.mod
cadecay.mod *
gabaa.mod
gabab.mod
HH2.mod *
Ih.mod *
IM.mod
IT.mod *
IT2.mod *
ITREcustom.mod
kleak.mod *
vecevent.mod
Fsinglecell.oc
Fspikewave.oc
membrane_potential_heat_plot.py
mosinit.hoc *
RE.tem
rundemo.hoc
screenshot1.png
screenshot2.png
screenshot3.png
sIN.tem
sPY.tem
TC.tem
                            
/*----------------------------------------------------------------------------

Implementation of:
Destexhe, Alain, "Spike-and-Wave Oscillations Based on the Properties of GABAB Receptors"
The Journal of Neuroscience, November 1, 1998, 18(21):9099–9111

Created by modifying the file Fspin.oc in Thalamocortical and Thalamic Reticular Network (Destexhe et al 1996, ModelDB Acession:3343) to create the four layer network described in the paper listed above, with various other functions added.

To vary cortical gabaa, change gabaapercent variable (ranges 0 - 100).  To vary T-type calcium channel parameters, change taubase_it2, shift_it2, or gcabar_it2, or similar variables for _itrecustom in the file RE.tem.

Andrew Knox 2014


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

ncorticalcells = 100
nthalamiccells = 100
narrowdiam = 5
widediam = 10

//change trans time if starting from a saved state (non-zero loads saved state)
trans = 0 // 10000
Dt = 0.1
dt = 0.1			// must be submultiple of Dt
npoints = 30000
stimtime = 10050			
randomstim = 0

showLFP = 0
fieldlower = 30
fieldupper = 70
fielddist = 50

watchneuron = 11 
axondelay = 0

gabaapercent = 0
gababpercent = 100

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

load_file("nrngui.hoc")

objectvar g[20]			// max 20 graphs
objectvar g1, g2, g3, g4
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(0)
	g[ii].view(0,1,0,1, int(ii/2)*550+80, ii%2*450+100, 400, 300)
	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 addfieldgraph() { local ii	// hack solution for a better positioned graph of electrical field
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph(0)
	g[ii].view(0,1,0,1, 650, 0, 800, 300)
	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.5,$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()		// create control menu
  ismenu=1
}

objref membranedatafile
membranedatafile = new File()



//----------------------------------------------------------------------------
// global simulation variables
//----------------------------------------------------------------------------

objectvar PY[ncorticalcells]		// create PY cells
objectvar PYVtrace[ncorticalcells]
objectvar IN[ncorticalcells]		// create IN cells
objectvar INVtrace[ncorticalcells]
objectvar RE[nthalamiccells]		// create TC cells
objectvar REVtrace[nthalamiccells]
objectvar TC[nthalamiccells]		// create RE cells
objectvar TCVtrace[nthalamiccells]

//stuff for raster plot
objref tobj, nil

objref PYtimevec, PYidvec, INtimevec, INidvec, REtimevec, REidvec, TCtimevec, TCidvec, recncs, recveclist

recncs = new List()
recveclist = new List()

field = 0

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

objectvar Sim			// create vector of simulation points
Sim = new Vector(npoints)

tstart = trans
tstop = trans + npoints * Dt
runStopAt = tstop
steps_per_ms = 1/Dt
celsius = 36
v_init = -70	

objectvar randvolt
randvolt = new Random()
randvolt.uniform(-80,-65)


objref lstate 
objref statefile

proc init () {
  finitialize(v_init)

//Assign random starting membrane voltages - doesn't really change anything if you let model come to steady state
//  for i=0,ncorticalcells-1 {
//    PY[i].soma.v(0.5) = randvolt.repick()
//    IN[i].soma.v(0.5) = randvolt.repick()
//  }
//  for i=0,nthalamiccells-1 {
//    TC[i].soma.v(0.5) = randvolt.repick()
//    RE[i].soma.v(0.5) = randvolt.repick()
//  }

  if (trans > 0)  {
    print "Setting trans to a non-zero value implies there is a state file to be loaded\n"
    lstate = new SaveState()
    statefile = new File()
    statefile.ropen("state_data.txt")
    lstate.fread(statefile)
    lstate.restore()
  }

  // params: RERE, RETCa,RETCb,TCRE, PYPY, PYIN, INPYa, INPYb, PYRE, PYTC, TCPY, TCIN
  assign_synapses(0.2,0.02,0.04,0.2,0.6,0.2,gabaapercent/100*0.15,gababpercent/100*0.03,1.2,0.01,1.2,0.4)  

  fcurrent()
  cvode.re_init()
  frecord_init()
}




//----------------------------------------------------------------------------
//  Create Cells
//----------------------------------------------------------------------------

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

load_file("TC.tem")		// read geometry file
for i=0,nthalamiccells-1 {
  TC[i] = new sTC()
}

load_file("RE.tem")		// read geometry file
for i=0,nthalamiccells-1 {
  RE[i] = new sRE()
}

load_file("sPY.tem")		// read geometry file
for i=0,ncorticalcells-1 {
  PY[i] = new sPY()
}

load_file("sIN.tem")		// read geometry file
for i=0,ncorticalcells-1 {
  IN[i] = new sIN()
}


//----------------------------------------------------------------------------
// set up recording stuff
//----------------------------------------------------------------------------

    PYtimevec = new Vector()
    PYidvec = new Vector()
    INtimevec = new Vector()
    INidvec = new Vector()
    REtimevec = new Vector()
    REidvec = new Vector()
    TCtimevec = new Vector()
    TCidvec = new Vector()

for i=0,ncorticalcells-1 {  
    PYVtrace[i] = new Vector()
    PYVtrace[i].record(&PY[i].soma.v(0.5))
    PY[i].soma tobj = new NetCon(&v(0.5), nil)
    tobj.record(PYtimevec, PYidvec, i+1) // so all the spike rasters lie above the x axis
    recncs.append(tobj)
    INVtrace[i] = new Vector()
    INVtrace[i].record(&IN[i].soma.v(0.5))
    IN[i].soma tobj = new NetCon(&v(0.5), nil)
    tobj.record(INtimevec, INidvec, i+1) // so all the spike rasters lie above the x axis
    recncs.append(tobj)
}


for i=0,nthalamiccells-1 {
    REVtrace[i] = new Vector()
    REVtrace[i].record(&RE[i].soma.v(0.5))
    RE[i].soma tobj = new NetCon(&v(0.5), nil)
    tobj.record(REtimevec, REidvec, i+1) // so all the spike rasters lie above the x axis
    recncs.append(tobj)
    TCVtrace[i] = new Vector()
    TCVtrace[i].record(&TC[i].soma.v(0.5))
    TC[i].soma tobj = new NetCon(&v(0.5), nil)
    tobj.record(TCtimevec, TCidvec, i+1) // so all the spike rasters lie above the x axis
    recncs.append(tobj)

}


//----------------------------------------------------------------------------
//   Set Up Synapses / Network Connectivity
//----------------------------------------------------------------------------

print " "
print "<<==================================>>"
print "<<     CREATE SYNAPTIC RECEPTORS    >>"
print "<<==================================>>"
print " "

func ncon() { local nc		// function to get the number of connections 
				// argument: interaction diameter
   nc = 2 * $1 + 1
   if(nc>ncorticalcells) nc = ncoricalcells
   return nc
}


//----------------------------------------------------
//  Glutamate AMPA receptors in synapses from TC to RE
//----------------------------------------------------

diamTCRE = narrowdiam		// diameter of connectivity for TC->RE
nTCRE = ncon(diamTCRE)	// nb of RE cells recampapost, PYlist, TClisteiving synapses from one TC cell

for i=0,nthalamiccells-1 {
   for j=i-diamTCRE,i+diamTCRE {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > nthalamiccells-1) jbound = 2 * nthalamiccells - jbound - 1

	// presynaptic is TC[i], postsynaptic is RE[j] 
	TC[i].soma RE[jbound].TClist.append(new NetCon(&v(0.5), RE[jbound].ampapost, 0, axondelay, 1))
   }
}

print " "
print "<< ",RE[0].TClist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM TC TO RE >>"
print " "


//---------------------------------------
//  GABAa receptors in intra-RE synapses
//---------------------------------------

diamRERE = narrowdiam		// diameter of connectivity for RE->RE
nRERE = ncon(diamRERE)	// nb of RE cells receiving synapses from one RE cell

for i=0,nthalamiccells-1 {
   for j=i-diamRERE,i+diamRERE {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > nthalamiccells-1) jbound = 2 * nthalamiccells - jbound - 1

	// presynaptic is RE[i], postsynaptic is RE[j] 
	RE[i].soma RE[jbound].REgabaalist.append(new NetCon(&v(0.5), RE[jbound].gabaapost, 0, axondelay, 1))
   }
}
	
print " "
print "<< ",RE[0].REgabaalist.count()," GABAa-MEDIATED SYNAPTIC CONTACTS FROM RE TO RE >>"
print " "


//--------------------------------------------------
//  GABAa receptors in synapses from RE to TC cells
//--------------------------------------------------

diamRETC = narrowdiam		// diameter of connectivity from RE->TC
nRETC = ncon(diamRETC)	// nb of RE cells receiving synapses from one TC cell

for i=0,nthalamiccells-1 {
   for j=i-diamRETC,i+diamRETC {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > nthalamiccells-1) jbound = 2 * nthalamiccells - jbound - 1

	// presynaptic is RE[i], postsynaptic is TC[j] 
	RE[i].soma TC[jbound].REgabaalist.append(new NetCon(&v(0.5), TC[jbound].gabaapost, 0, axondelay, 1))
   }
}
	
print " "
print "<< ",TC[jbound].REgabaalist.count()," GABAa-MEDIATED SYNAPTIC CONTACTS FROM RE TO TC >>"
print " "



//--------------------------------------------------
//  GABAb receptors in synapses from RE to TC cells
//--------------------------------------------------

// use same diameters and connectivity as GABAa receptors (colocalized)
objectvar gababsyn

for i=0,nthalamiccells-1 {
   for j=i-diamRETC,i+diamRETC {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > nthalamiccells-1) jbound = 2 * nthalamiccells - jbound - 1

	// presynaptic is RE[i], postsynaptic is TC[j]
	// ***Note: GABAb synapses are implemented as a list of individual synapses (in contrast to other synapse types), and so are created here
	gababsyn = new GABAb_S()
 	TC[jbound].soma gababsyn.loc(0.5)
	TC[jbound].gababpost.append(gababsyn)
	RE[i].soma TC[jbound].REgabablist.append(new NetCon(&v(0.5), gababsyn, 0, axondelay, 1))	
   }
}
	
print " "
print "<< ",TC[0].REgabablist.count()," GABAb-MEDIATED SYNAPTIC CONTACTS FROM RE TO TC >>"
print " "


//-----------------------------------------------------
//  Glutamate AMPA receptors in synapses from PY to PY
//-----------------------------------------------------

diamPYPY = narrowdiam		// diameter of connectivity for PY->PY
nPYPY = ncon(diamTCRE)	// nb of PY cells receiving synapses from one PY cell

for i=0,ncorticalcells-1 {
   for j=i-diamPYPY,i+diamPYPY {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > ncorticalcells-1) jbound = 2 * ncorticalcells - jbound - 1
	
	// presynaptic is PY[i], postsynaptic is PY[j] 
	//if (i != jbound) PY[i].soma PY[jbound].PYlist.append(new NetCon(&v(0.5), PY[jbound].ampapostPY, 0, axondelay, 1))  //using this line rather than the line below eliminates self connections
	PY[i].soma PY[jbound].PYlist.append(new NetCon(&v(0.5), PY[jbound].ampapostPY, 0, axondelay, 1))	
   }
}

print " "
print "<< ",PY[0].PYlist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM PY TO PY >>"
print " "



//------------------------------------------------------
//  Glutamate AMPA receptors in synapses from PY to IN
//------------------------------------------------------

diamPYIN = narrowdiam		// diameter of connectivity for PY->IN
nPYIN = ncon(diamPYIN)	

for i=0,ncorticalcells-1 {
   for j=i-diamPYIN,i+diamPYIN {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > ncorticalcells-1) jbound = 2 * ncorticalcells - jbound - 1

	// presynaptic is PY[i], postsynaptic is IN[j] 
	PY[i].soma IN[jbound].PYlist.append(new NetCon(&v(0.5), IN[jbound].ampapost, 0, axondelay, 1))	
   }
}

print " "
print "<< ",IN[0].PYlist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM PY TO IN >>"
print " "



//--------------------------------------------------
//  GABAa receptors in synapses from IN to PY cells
//--------------------------------------------------

diamINPY = narrowdiam		// diameter of connectivity from IN->PY
nINPY = ncon(diamINPY)	

for i=0,ncorticalcells-1 {
   for j=i-diamINPY,i+diamINPY {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > ncorticalcells-1) jbound = 2 * ncorticalcells - jbound - 1

	// presynaptic is IN[i], postsynaptic is PY[j] 
	IN[i].soma PY[jbound].INgabaalist.append(new NetCon(&v(0.5), PY[jbound].gabaapost, 0, axondelay, 1))	
   }
}
	
print " "
print "<< ",PY[0].INgabaalist.count()," GABAa-MEDIATED SYNAPTIC CONTACTS FROM IN TO PY >>"
print " "


//--------------------------------------------------
//  GABAb receptors in synapses from IN to PY cells
//--------------------------------------------------

// use same diameters and connectivity as GABAa receptors (colocalized)

for i=0,ncorticalcells-1 {
   for j = i-diamINPY, i+diamINPY {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > ncorticalcells-1) jbound = 2 * ncorticalcells - jbound - 1

	// presynaptic is IN[i], postsynaptic is PY[j]
	// ***Note: GABAb synapses are implemented as a list of individual synapses (in contrast to other synapse types), and so are created here
	gababsyn = new GABAb_S()
 	PY[jbound].soma gababsyn.loc(0.5)
	PY[jbound].gababpost.append(gababsyn) 
	IN[i].soma PY[jbound].INgabablist.append(new NetCon(&v(0.5), gababsyn, 0, axondelay, 1))	
   }
}
	
print " "
print "<< ",PY[0].INgabablist.count()," GABAb-MEDIATED SYNAPTIC CONTACTS FROM IN TO PY >>"
print " "



//----------------------------------------------------------------------------
//  Glutamate AMPA receptors in synapses from PY to RE
//----------------------------------------------------------------------------

diamPYRE = widediam		
nPYRE = ncon(diamPYRE)	
divergence = ncorticalcells / nthalamiccells

for (i=0;i<=ncorticalcells-1;i+=divergence) {
   for j=i/divergence-diamPYRE,i/divergence+diamPYRE {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > nthalamiccells-1) jbound = 2 * nthalamiccells - jbound - 1

	// presynaptic is PY[i], postsynaptic is RE[j] 
	PY[i+divergence/2].soma RE[jbound].PYlist.append(new NetCon(&v(0.5), RE[jbound].ampapost, 0, axondelay, 1))	
   }
}

print " "
print "<< ",RE[0].PYlist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM PY TO RE >>"
print " "


//----------------------------------------------------------------------------
//  Glutamate AMPA receptors in synapses from PY to TC
//----------------------------------------------------------------------------

diamPYTC = widediam	
nPYTC = ncon(diamPYTC)	

for (i=0;i<=ncorticalcells-1;i+=divergence) {
   for j=i/divergence-diamPYTC,i/divergence+diamPYTC {
	jbound = j
	if (jbound < 0) jbound = abs(j) - 1
	if (jbound > nthalamiccells-1) jbound = 2 * nthalamiccells - jbound - 1

	// presynaptic is PY[i], postsynaptic is RE[j] 
	PY[i+divergence/2].soma TC[jbound].PYlist.append(new NetCon(&v(0.5), TC[jbound].ampapost, 0, axondelay, 1))	
   }
}

print " "
print "<< ",TC[0].PYlist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM PY TO TC >>"
print " "

	


//----------------------------------------------------------------------------
//  Glutamate AMPA receptors in synapses from TC to PY
//----------------------------------------------------------------------------

diamTCPY = widediam		
nTCPY = ncon(diamTCPY)	

for i=0,nthalamiccells-1 {
   for (j=divergence*(i-diamTCPY);j<=divergence*(i+diamTCPY);j+=divergence) {
	jbound = j
	if (jbound < 0) jbound = abs(j) - divergence
	if (jbound > ncorticalcells-1) jbound = 2 * ncorticalcells - jbound - divergence

	// presynaptic is TC[i], postsynaptic is PY[j] 
	TC[i].soma PY[jbound+divergence/2].TClist.append(new NetCon(&v(0.5), PY[jbound+divergence/2].ampapostTC, 0, axondelay, 1))
   }
}

print " "
print "<< ",PY[0].TClist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM TC TO RE >>"
print " "



//----------------------------------------------------------------------------
//  Glutamate AMPA receptors in synapses from TC to IN
//----------------------------------------------------------------------------

diamTCIN = widediam		// diameter of connectivity for TC->IN
nTCIN = ncon(diamTCIN)	

for i=0,nthalamiccells-1 {
   for (j=divergence*(i-diamTCIN);j<=divergence*(i+diamTCIN);j+=divergence) {
	jbound = j
	if (jbound < 0) jbound = abs(j) - divergence
	if (jbound > ncorticalcells-1) jbound = 2 * ncorticalcells - jbound - divergence

	// presynaptic is TC[i], postsynaptic is IN[j] 
	TC[i].soma IN[jbound+divergence/2].TClist.append(new NetCon(&v(0.5), IN[jbound+divergence/2].ampapost, 0, axondelay, 1))
   }
}

print " "
print "<< ",IN[0].TClist.count()," AMPA-MEDIATED SYNAPTIC CONTACTS FROM TC TO IN >>"
print " "

	
//Synapse parameters
	Alpha_AMPA_S = 0.94		// kinetics from simplex with short pulses
	Beta_AMPA_S = 0.18 //0.18
	Cmax_AMPA_S = 0.5
	Cdur_AMPA_S = 0.3
	Erev_AMPA_S = 0

	Alpha_GABAa_S = 20		// from diffusion model
	Beta_GABAa_S = 0.162
	Cmax_GABAa_S = 0.5		// short pulses
	Cdur_GABAa_S = 0.3
	Erev_GABAa_S = -85		// Rinzel's Erev

	K1_GABAb_S	= 0.09		//	(/ms mM) forward binding to receptor
	K2_GABAb_S	= 0.0012	//	(/ms)	backward (unbinding)of receptor
	K3_GABAb_S	= 0.18 		//	(/ms)	rate of G-protein production
	K4_GABAb_S	= 0.034		//	(/ms)	rate of G-protein decay  -  larger number = slower decay?
	KD_GABAb_S	= 100		//	dissociation constant of K+ channel
	n_GABAb_S	= 4		//	nb of binding sites of G-protein on K+
	Erev_GABAb_S	= -95		//	(mV)	reversal potential (E_K)
	Cmax_GABAb_S = 0.5		// short pulses
	Cdur_GABAb_S = 0.3



proc assign_synapses() {	// procedure to assign syn conductances 
				// params: 1=intraRE (RERE), 2=GABAa in TC (RETC),
				// 3=GABAb in TC, 4=AMPA in RE (TCRE)
				// 5=PYPY 6=PYIN 7=GABAa INPY
				// 8=GABAb INPY 9=PYRE 10=PYTC
				// 11=TCPY 12=TCIN
  print "nRERE:", nRERE, "nRETC:", nRETC, "nTCRE:", nTCRE, "nPYPY:", nPYPY, "nINPY:", nINPY, "nPYRE:", nPYRE, "nPYTC:", nPYTC, "nTCPY:", nTCPY, "nTCIN:", nTCIN
  for i=0, nthalamiccells-1 {
    for j=0, RE[i].REgabaalist.count()-1 {
	RE[i].REgabaalist.object(j).weight = $1 / nRERE               
    }
    for j=0, TC[i].REgabaalist.count()-1 {
	TC[i].REgabaalist.object(j).weight = $2 / nRETC
    }
    for j=0, TC[i].gababpost.count()-1 {
	TC[i].gababpost.object(j).gmax = $3 / nRETC
	//TC[i].REgabablist.object(j).weight = $3 / nRETC
    }
    for j=0, RE[i].TClist.count()-1 {
	RE[i].TClist.object(j).weight = $4 / nTCRE
    }
 
    for j=0, RE[i].PYlist.count()-1 {  		
	RE[i].PYlist.object(j).weight = $9 / nPYRE
    }
    for j=0, TC[i].PYlist.count()-1 {
	TC[i].PYlist.object(j).weight = $10 / nPYTC
    }
  }
  for i=0, ncorticalcells-1 {
    for j=0, PY[i].PYlist.count()-1 {
	PY[i].PYlist.object(j).weight = $5 / nPYPY
	print "PY", i, j, PY[i].PYlist.object(j).weight
    }
    for j=0, IN[i].PYlist.count()-1 {
	IN[i].PYlist.object(j).weight = $6 / nPYIN
    }
    for j=0, PY[i].INgabaalist.count()-1 {
	PY[i].INgabaalist.object(j).weight = $7 / nINPY
    }
    for j=0, PY[i].gababpost.count()-1 {
	PY[i].gababpost.object(j).gmax = $8 / nINPY
	//PY[i].INgabablist.object(j).weight = $8 / nINPY
    }
    for j=0, PY[i].TClist.count()-1 {			
	PY[i].TClist.object(j).weight  = $11 / nTCPY
    }
    for j=0, IN[i].TClist.count()-1 {			
	IN[i].TClist.object(j).weight = $12 / nTCIN
    }
  }
}


//---------------------------------------------------------------------------
//   Assign potassium leak and Ih current strengths in TC cells 
//---------------------------------------------------------------------------

objectvar rgh, rk1
rgh = new Random()
rk1 = new Random()
rgh.normal(17.5,0.0008)        //random number generator behaves weirdly for very small numbers, so multiply by 10^-6 below
rk1.normal(40,0.003)

// setup TC cells in resting mode (no spontaneous oscillation)
for i=0,nthalamiccells-1 { 
  TC[i].soma.ghbar_iar = rgh.repick() * 10^-6   
  TC[i].kl.gmax = rk1.repick() * 10^-4         
  print "TC(",i,") gh:",TC[i].soma.ghbar_iar," gmax:",TC[i].kl.gmax

}



//----------------------------------------------------------------------------
//  set up current stimulus to cells
//----------------------------------------------------------------------------

smallPY = 1
mediumPY = 0
largePY = 0
largePYhole = 0
smallPYhole = 0
mediumPYoffset = 7
largePYoffset = 2

objectvar PYstim[ncorticalcells/5]		
objectvar jitter

if (largePY==1) {
  jitter = new Random()
  jitter.uniform(-1,1)

  for i = 0, (ncorticalcells/5)-1 PY[i*(ncorticalcells/20)+largePYoffset].soma { 
	PYstim[i] = new IClamp(0.5)
	PYstim[i].dur = 100
	PYstim[i].amp = 0.7

	if(randomstim==1) {
	  PYstim[i].del = stimtime + jitter.repick()
	} else {
	  PYstim[i].del = stimtime
	}
  }
  if (largePYhole==1)  {
    PYstim[8].amp = 0
    PYstim[9].amp = 0
    PYstim[11].amp = 0
    PYstim[12].amp = 0
  }
  if (smallPYhole==1) {
    PYstim[9].amp = 0
    PYstim[10].amp = 0
  }
}

if (smallPY==1) {
  //for i = 0, (ncorticalcells/20)-1 PY[i*(ncorticalcells/5)+9].soma {
  for i = 0, (ncorticalcells/20)-1 PY[i*(ncorticalcells/5-1)+11].soma {
	PYstim[i] = new IClamp(0.5)
	PYstim[i].del = stimtime
	PYstim[i].dur = 100
	PYstim[i].amp = 0.7
  }
}

if (mediumPY==1) {
  for i = 0, (ncorticalcells/10)-3 PY[i*(ncorticalcells/10+2)+mediumPYoffset].soma { //most were done with +1
	PYstim[i] = new IClamp(0.5)
	PYstim[i].del = stimtime
	PYstim[i].dur = 100
	PYstim[i].amp = 0.7
  }
  if (smallPYhole==1) {
    PYstim[5].amp=0
  }
}

/*
//use this code to give stim to RE neurons
objectvar REstim[nthalamiccells/5]
for i = 0, (nthalamiccells/5)-1 RE[i*(nthalamiccells/20)+nthalamiccells/10].soma {
	REstim[i] = new IClamp(0.5)
	REstim[i].del = stimtime
	REstim[i].dur = 100
	REstim[i].amp = 0
}
*/

/*
//use this code to give stim to TC neurons
objectvar TCstim[nthalamiccells/5]
//for i = 0, (nthalamiccells/20)-1 TC[i*(nthalamiccells/5)+nthalamiccells/10].soma {
//for i = 0, (nthalamiccells/10)-1 TC[i*(nthalamiccells/10)+5].soma {
for i = 0, (nthalamiccells/5)-1 TC[i*(nthalamiccells/20-1)+10].soma {
//for i = 0, 9 TC[i*9+9].soma {
	TCstim[i] = new IClamp(0.5)
	TCstim[i].del = stimtime
	TCstim[i].dur = 100
	TCstim[i].amp = -0.1
}
*/




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

//addgraph("tcB[0][0].g",0,0.05)
//addgraph("TC[0].soma.o2_iar",0,1)
//addgraph("TC[0].soma.p1_iar",0,1)
//addgraph("Vampa[50]",-70,15)
//addgraph("TC[50].soma.ina",-.02,.001)

//addgraph("TC[watchneuron].soma.ica",-0.01,0)
//addgraph("TC[watchneuron].soma.ih",-.001,.001)
//addgraph("PY[watchneuron].ampapostPY.i", -20, 0)
//addgraph("PY[watchneuron].ampapost.synon", 0, 2)
//addgraph("PY[watchneuron].soma.ma_hh2",-.002,1)
//addgraph("PY[watchneuron].soma.mb_hh2",-.002,1)
//addgraph("PY[watchneuron].soma.ik_im",-.002,.1)
//addgraph("TC[50].soma.ica",-.02,.001)

//addgraph("RE[50].soma.ica",-.02,.001)
//addgraph("RE[50].soma.ica_itrecustom",-.02,.001)


strdef gtxt
//for i=0,ncells-1 {
//for i=50,51 {
i=watchneuron
	sprint(gtxt,"TC[%d].soma.v(0.5)",i)
	addgraph(gtxt,-120,40)
	sprint(gtxt,"RE[%d].soma.v(0.5)",i)
	addgraph(gtxt,-120,40)
        sprint(gtxt,"PY[%d].soma.v(0.5)",i)
	addgraph(gtxt,-120,40)
	sprint(gtxt,"IN[%d].soma.v(0.5)",i)
	addgraph(gtxt,-120,40)
//}

if (showLFP) addfieldgraph("field",-5,5)


//-----------------------------------------------------------------------------
// Make raster plots
//-----------------------------------------------------------------------------


proc rasterplot() {
  print "plot:", tstart, tstop
  //PYtimevec.printf()
  //PYidvec.printf()
  plotlen = 3000 
  g1 = new Graph(0)
  g1.label(0.5, 0.95, "PY")
  g1.view(tstart,0,plotlen,100,100,500,400,300)
  PYidvec.mark(g1, PYtimevec,"o",2)
  g2 = new Graph(0)
  g2.label("IN")
  g2.view(tstart,0,plotlen,100,650,500,400,300)
  INidvec.mark(g2, INtimevec,"o",2)
  g3 = new Graph(0)
  g3.label(0.5, 0.95,"RE")
  g3.view(tstart,0,plotlen,100,650,500,400,300)
  REidvec.mark(g3, REtimevec,"o",2)
  g4 = new Graph(0)
  g4.label(0.5, 0.95,"TC")
  g4.view(tstart,0,plotlen,100,1200,500,400,300)
  TCidvec.mark(g4, TCtimevec,"o",2)
}



//-----------------------------------------------------------------------------
// Write recorded membrane data to file
//-----------------------------------------------------------------------------

proc writedatafile()  {

  membranedatafile.wopen("membrane_data.txt")


//only take every other time point to avoid annoying memory errors.
  membranedatafile.printf("%d\n",TCVtrace[0].size())
  membranedatafile.printf("%d\n",ncorticalcells)
  for i=0,nthalamiccells-1 {
    TCVtrace[i].printf(membranedatafile)
  }
  for i=0,nthalamiccells-1{
    REVtrace[i].printf(membranedatafile)
  }
  for i=0,ncorticalcells-1{
    PYVtrace[i].printf(membranedatafile) 
  }
  for i=0,ncorticalcells-1{
    INVtrace[i].printf(membranedatafile)
  }

  membranedatafile.close()
}

//-------------------------------------------------------------------------------
// Save state - state will automatically load if trans is set to a non-zero value
//-------------------------------------------------------------------------------

objref sstate, statefile 

proc writestate()  {
   sstate = new SaveState()
   sstate.save()
   statefile = new File()
   statefile.wopen("state_data.txt")
   sstate.fwrite(statefile)
}


//---------------------------------------------------------------------------
//  Code for dealing with field potentials
//---------------------------------------------------------------------------

func xfield()  {local i, j, tmp, Ni, Re, x, cw   
//$1 is Re, $2 is distance from nearest neuron, $3 is neuron to which probe is closest
//assumes field comes from PY neurons, which are in a row spaced 20um apart
   total_field = 0
   Re = $1
   x = $2
   Ni = $3

   for i = fieldlower, fieldupper {
	tmp = 0
	for j = 0, PY[i].gababpost.count()-1 {
	   tmp += PY[i].gababpost.object(j).i
	}

	tmp += PY[i].gabaapost.i  + PY[i].soma.ik_im + PY[i].ampapostPY.i + PY[i].ampapostTC.i
	total_field += tmp / sqrt(x^2 + 400*(Ni-i)^2)
   }
   return total_field * Re/4/PI
}

proc advance() {
   fadvance()
   field = xfield(230,fielddist,watchneuron)
}


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

access TC[0].soma

proc text() {
  sprint(gtxt,"%d RE and %d TC cell",nthalamiccells,nthalamiccells)
  addtext(gtxt)
  sprint(gtxt,"Membrane level: kleak=%g",TC.kl.gmax)
  addline(gtxt)
  sprint(gtxt,"Ih: g=%g, ginc=%g, nca=%g, k2=%g, cac=%g", \
  ghbar_iar,ginc_iar,nca_iar,k2_iar,cac_iar)
  addline(gtxt)
  sprint(gtxt,"Ih: nexp=%g, Pc=%g, k4=%g, taum=%g", \
  nexp_iar,Pc_iar,k4_iar,taum_iar)
  addline(gtxt)
  sprint(gtxt,"GABAa_S: Alpha=%g, Beta=%g, Cdur=%g, Cmax=%g", \
  Alpha_GABAa_S,Beta_GABAa_S,Cmax_GABAa_S,Cdur_GABAa_S)
  addline(gtxt)
  sprint(gtxt,"GABAb_S: K1=%g, K2=%g, K3=%g, K4=%g, KD=%g", \
  K1_GABAb_S,K2_GABAb_S,K3_GABAb_S,K4_GABAb_S,KD_GABAb_S)
  addline(gtxt)
}

print " "
print "Use procedure text() to create a new window with actual parameters"
print "Use procedure assign_synapses() to change synaptic conductances"
print " "

Loading data, please wait...