// 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 9F
// in Ovsepian et al. (2013)


strdef strFilePrefix

strFilePrefix = "Kdr60" // will be overwritten

load_file("nrngui.hoc")
load_file("DCN_params_axis.hoc")
load_file("DCN_morph.hoc")
load_file("DCN_mechs1.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 (ccip = 0; ccip < 2; ccip += 1) {

	SOMACIP1AMP  = -0.0947 //-0.02*ccip - 0.200 // nA

	sprint(strFilePrefix, "SomaCip_m%d_e-1pA_Kdrblock%d", int(-10000*SOMACIP1AMP), int(100*Kdrblock))
	print "writing output to file with prefix ", strFilePrefix
  
	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()