/*----------------------------------------------------------------------------
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 " "
