// CN model used in Saak V Ovsepian, Volker Steuber, Marie Le
// Berre, Liam O'Hara, Valerie B O'Leary, and J. Oliver Dolly
// (2013). A Defined Heteromeric KV1 Channel Stabilizes the
// Intrinsic Pacemaking and Regulates the Efferent Code of Deep
// Cerebellar Nuclear Neurons to Thalamic Targets. Journal of
// Physiology (epub ahead of print).
//
// written by Johannes Luthman, modified by Volker Steuber
//
// main simulation script that replicates Figure 9A,B
// in Ovsepian et al. (2013)
strdef strFilePrefix
Kdrblock = 0.6 // will be overwritten
strFilePrefix = "Kdr60"
load_file("nrngui.hoc")
load_file("DCN_params.hoc")
load_file("DCN_morph.hoc")
load_file("DCN_mechs2.hoc")
objref oRndInh, oRndExc
objref gammaStimPC[PCtoDCNconvergence]
objref netconPC[INHTOTALSYNAPSES]
// Declare instances of the GammaStim objects for excitatory synapses
// (1 GammaStim activates 1 AMPA + 1 fNMDA + 1 sNMDA) and the corresponding
// NetCon objects (=1 each for AMPA, fNMDA, and sNMDA).
objref gammaStimExc[EXCTOTALSYNAPSES]
objref netconExc[3 * EXCTOTALSYNAPSES]
///////////////////////////////////////////////////////////////
proc DCNloop() {
for (cKdrblock = 0; cKdrblock < 9; cKdrblock += 1) {
Kdrblock = 0.6 + cKdrblock*0.1
sprint(strFilePrefix, "I020pA_Kdr%d", int(100*Kdrblock))
print "writing output to file with prefix ", strFilePrefix
gfKdrsoma = 1.5e-2*QdTconductances*Kdrblock
gfKdraxHill = 2*gfKdrsoma*Kdrblock
gfKdraxIniSeg = 2*gfKdrsoma*Kdrblock
gfKdrpDend = 0.6*gfKdrsoma*Kdrblock
gsKdrsoma = 1.25e-2*QdTconductances*Kdrblock
gsKdraxHill = 2*gsKdrsoma*Kdrblock
gsKdraxIniSeg = 2*gsKdrsoma*Kdrblock
gsKdrpDend = 0.6*gsKdrsoma*Kdrblock
DCNmechs()
load_file("DCN_recording.hoc")
setupSimulation()
runSimulation()
}
quit()
}
proc DCNrun() {
load_file("DCN_recording.hoc")
setupSimulation()
runSimulation()
quit()
}
proc setupSimulation() {
// Work around the error that sets in when both tTraceStart and tTraceStop = 0
// (the latter setting is used to not record any traces).
if (tTraceStop[0] == 0) {
tTraceStart[0] = 1
}
// Set up size of recording vectors.
nRecVectElements = int(runTime/recInterval) + 2 //+1 is ok for recInterval=100us,
// but +2 is needed for recInterval=50us.
stepsPerRec = int(recInterval/dt)
// Set up a different size for the trace vectors.
if (tTraceStop[0] > 0) {
recordTraces = 1
sizeVectorOfTrace = getSizeOfTraceVectors()
} else {
recordTraces = 0
}
} // end of proc setupSimulation()
proc runSimulation() {
// Set up the excitatory synapses
oRndExc = new Random()
oRndExc.ACG(randomiserSeed)
// Add GammaStims and connect them with the synaptic mechanisms using NetCons.
// As an ARTIFICIAL_CELL, it doesn't matter where a GammaStim.mod is placed.
// For convenience, I place all the GammaStims in the soma.
for (c=0; c < EXCTOTALSYNAPSES; c+=1) {
soma gammaStimExc[c] = new GammaStim(0.5)
gammaStimExc[c].start = 0
gammaStimExc[c].noise = noiseFractionExcSyn
gammaStimExc[c].duration = runTime
gammaStimExc[c].order = gammaOrderExc
gammaStimExc[c].refractoryPeriod = refractoryPeriodExc
}
// Set up the excitatory NetCons.
ncIndex = 0
for (c=0; c < EXCTOTALSYNAPSES; c+=1) {
// Create three NetCons for each excitatory GammaStim to connect it
// to the ampa and the two nmda receptors.
netconExc[ncIndex] = new NetCon(gammaStimExc[c], ampa[c])
netconExc[ncIndex].threshold = 0 //mV
netconExc[ncIndex].weight = gAMPA
netconExc[ncIndex+1] = new NetCon(gammaStimExc[c], fnmda[c])
netconExc[ncIndex+1].threshold = 0 //mV
netconExc[ncIndex+1].weight = gfNMDA
netconExc[ncIndex+2] = new NetCon(gammaStimExc[c], snmda[c])
netconExc[ncIndex+2].threshold = 0 //mV
netconExc[ncIndex+2].weight = gsNMDA
ncIndex = ncIndex + 3
}
// Set up the inhibitory synapses
// For synapses not incorporating synaptic depression, set the max GABA
// conductance to equal what is reached for depressed synapses on
// completely regular input at the used input frequency.
if (useGABAsyndep == 0) {
// The following calculation is the same as in DCN_mechs.hoc where it's explained:
initDeprLevel = 0.08 + 0.60*exp(-2.84*inhibitoryHz) + 0.32*exp(-0.02*inhibitoryHz)
gGABA = gGABA*initDeprLevel
}
oRndInh = new Random()
oRndInh.ACG(randomiserSeed)
for (c=0; c<PCtoDCNconvergence; c+=1) {
soma gammaStimPC[c] = new GammaStim(0.5)
gammaStimPC[c].start = 0
gammaStimPC[c].noise = noiseFractionInhSyn
gammaStimPC[c].duration = runTime
gammaStimPC[c].order = gammaOrderPC
gammaStimPC[c].refractoryPeriod = refractoryPeriodPC
}
// Set up the GABA NetCons.
gsIndex = 0
counterOfNetCons = 0
for (cGABA=0; cGABA < INHTOTALSYNAPSES; cGABA=cGABA+1) {
netconPC[cGABA] = new NetCon(gammaStimPC[gsIndex], gaba[cGABA])
netconPC[cGABA].threshold = 0 //mV
netconPC[cGABA].weight = gGABA
// If we've used all synapses that fit onto one GammaStim, then start
// using the next GammaStim.
if (counterOfNetCons == (nDCNsynsPerPC-1)) {
gsIndex += 1
counterOfNetCons = 0
} else {
counterOfNetCons+=1
}
}
// All the artifical cell elements have been set up. Seed them. Calling
// .seed() affects the event streams generated by all GammaStims,
// and thus needs to be done for only one of them.
gammaStimExc[0].seed(randomiserSeed)
iRecTimeVolt = 0 //index of spike recording vectors (time and volt)
// If the user has set to record volt/current traces then setup index for that.
if (tTraceStop[0] > 0) {
iRecTrace = 0
}
// Counter for when to record variables
ndtSinceRecording = stepsPerRec
InstantiateRecObjects() // procedure in file DCN_recording.hoc
SetupOutputFiles() // procedure in file DCN_recording.hoc
t = 0
print "initialising simulation"
finitialize (vInit)
// Set properties of the inhibitory GammaStim elements
if (inhibitoryHz > 0) {
intervalInh = 1000 / inhibitoryHz
delayInh = oRndInh.uniform(0, intervalInh)
} else {
// No inhibitory input is to be provided - set delay so the
// first spike occurs outside of the simulation time.
intervalInh = 1e5
delayInh = runTime
}
for (c=0; c < PCtoDCNconvergence; c+=1) {
gammaStimPC[c].interval = intervalInh
}
// Set delays of the inhibitory netCons.
counterOfNetCons = 0
for (cGABA=0; cGABA < INHTOTALSYNAPSES; cGABA=cGABA+1) {
netconPC[cGABA].delay = delayInh
if (counterOfNetCons == (nDCNsynsPerPC-1)) {
if (inhibitoryHz > 0) {
delayInh = oRndInh.repick()
}
counterOfNetCons = 0
} else {
counterOfNetCons+=1
}
}
// Excitatory synaptic inputs
if (excitatoryHz > 0) {
intervalExc = 1000 / excitatoryHz
delayExc = oRndExc.uniform(0, intervalExc)
} else {
// No excitatory input is to be provided - set delay so the
// first spike occurs outside of the simulation time.
intervalExc = 1e5
delayExc = runTime
}
for (c=0; c < EXCTOTALSYNAPSES; c+=1) {
gammaStimExc[c].interval = intervalExc
}
// Set delay for this stage.
for (c=0; c < (3 * EXCTOTALSYNAPSES); c=c+3) {
if (excitatoryHz > 0) {
delayExc = oRndExc.repick()
}
netconExc[c].delay = delayExc
netconExc[c+1].delay = delayExc
netconExc[c+2].delay = delayExc
}
// Step through the simulation
print "starting simulation"
while (t < runTime) {
// Write variables to vectors if the time interval since the last recording
// has been reached
if (ndtSinceRecording == stepsPerRec) {
writeToTimeAndVoltVectors() // procedure in file DCN_recording.hoc
iRecTimeVolt+=1
if (recordTraces == 1) {
if (isItTimeToSaveTrace() == 1) {
writeToTraceVectors()
iRecTrace = iRecTrace+1
}
}
ndtSinceRecording = 0
}
fadvance()
//print t
ndtSinceRecording+=1
}
print "done"
// Save any recordings of traces to file (they're usually set to only a small interval
// of the whole simulation, making it more efficient to save all of it here at the end.
if (recordTraces == 1) {
writeTracesToFile() // procedure in file DCN_recording.hoc
}
writeSpikeTimesToFile()
} // end of "proc runSimulation()".
func isItTimeToSaveTrace() { local subBoolean
// Check if the current time is within the time frame set in DCN_simulation.hoc
// for recording of the traces (there's no "exit for" in hoc..)
subBoolean = 0
for (c=0; c < nStepsSaveTrace; c+=1) {
if (t>=tTraceStart[c] && t<=tTraceStop[c]) {
subBoolean = 1
}
}
return subBoolean
}
func getSizeOfTraceVectors() { local subC, subSumIndeces
// For each step of recording the full voltage/current trace/s, calculate
// how many vector indeces will be needed for it. In the following,
// +1 is to get the extra index a exactly time
// tTraceStart[subC].
subSumIndeces = 0
for (subC=0; subC < nStepsSaveTrace; subC=subC+1) {
subSumIndeces = subSumIndeces + int((tTraceStop[subC] - tTraceStart[subC]) \
/ recInterval) + 2
}
return int(subSumIndeces)
}
func getPCtoDCNdelay() { local subTemp, subMean, subSD
subMean = $1
subSD = $2
if (subSD >= 0.000001) {
subTemp = oRndInh.repick()
while (subTemp < 0) {
subTemp = oRndInh.repick()
}
} else {
subTemp = subMean
}
return subTemp
}
//xpanel("Luthman et al. 2011")
// xbutton("DCNrun()","DCNrun()")
//xpanel()
//DCNrun()
DCNloop()
| 
