/* Author: Xu Zhang @UConn, Jan., 2020
This script implements real-time phase-locked tACS for essential tremor.
It must be run after INIT_ET_tACS, which generates the state file "currentState.dat"
and initializes recording data files in the "temporary" and "recordings_full" folders.
Each time this script is executed, it loads te state file "currentState.dat",
runs the simulation for one time step, then rewrites state and data.
*/

load_file("nrngui.hoc")

load_file("rngRead.hoc") // Load rng seeds
noiseSwitch=1 // Membrane and synaptic noises, 1/0 = turn on/off

// Loading individual cells
load_file("cell_ION.hoc")
load_file("cell_PC.hoc")
load_file("cell_DCN.hoc")
load_file("cell_TC.hoc")
load_file("cell_MC.hoc")

// Loading synaptic connections
load_file("syn_ION_PC.hoc")
load_file("syn_PC_DCN.hoc")
load_file("syn_DCN_ION.hoc")
load_file("syn_DCN_TC.hoc")
load_file("syn_TC_MC.hoc")
load_file("syn_MC_GrL.hoc")
load_file("syn_GrL_PC.hoc")

// One typical ET condition (tau = 12, R = 0.7)
PCDCNtaurnd = 682
PCDCNgrnd = 198
objref rngPCDCNtau, PCDCNg
rngPCDCNtau = new Random(PCDCNtaurnd)
rngPCDCNtau.normal(12,4)
PCDCNg = new Random(PCDCNgrnd)
PCDCNg.normal(0.7,0.0196)
for i = 0,199 {
	PC2DCN_syn[i].tau = rngPCDCNtau.repick()
	PC2DCN_syn[i].g = PC2DCN_syn[i].g*PCDCNg.repick()
}

// Push ION into oscillations; trigger tremor onset
objref stimION[40]
for i = 0,39 {
	IONcell[i] stimION[i] = new IClamp(0.5)
	stimION[i].amp = 1e-2
	stimION[i].del = 1000
	stimION[i].dur = 20
access DCNcell
}


objref f_trialname
f_trialname = new File()
f_trialname.ropen("trialname.txt")

strdef f_PCap, f_ION
strdef f_DCN1, f_DCN2, f_DCN3, f_DCN4, f_DCN5
strdef f_PCv
strdef f_Vim1, f_Vim2, f_Vim3, f_Vim4, f_Vim5
strdef f_sin
f_trialname.scanstr(f_PCap)
f_trialname.scanstr(f_ION)
f_trialname.scanstr(f_DCN1)
f_trialname.scanstr(f_DCN2)
f_trialname.scanstr(f_DCN3)
f_trialname.scanstr(f_DCN4)
f_trialname.scanstr(f_DCN5)
f_trialname.scanstr(f_PCv)
f_trialname.scanstr(f_Vim1)
f_trialname.scanstr(f_Vim2)
f_trialname.scanstr(f_Vim3)
f_trialname.scanstr(f_Vim4)
f_trialname.scanstr(f_Vim5)
f_trialname.scanstr(f_sin)
f_trialname.close()

objref f_sinparams
f_sinparams = new File()
f_sinparams.ropen("sinparams.txt")
ampparam = f_sinparams.scanvar()
phaseparam = f_sinparams.scanvar()
f_sinparams.close()

objref sinuStim_group[200]
for i = 0,199 {
PCcell[i] sinuStim_group[i] = new SinClamp(0.3)
sinuStim_group[i].phase = phaseparam
sinuStim_group[i].pkamp = ampparam
sinuStim_group[i].freq = 7
sinuStim_group[i].del = 0
sinuStim_group[i].dur = 1e5
}


v_init = -57
dt = 0.0125
tstop = 0.25
init()

objref ss,f_ss
ss = new SaveState()
f_ss = new File("currentState.dat")
ss.fread(f_ss)
ss.restore()
f_ss.close()

for i = 0,199 {
sinuStim_group[i].phase = phaseparam // 7/4*PI
}

// Record PC
objref apc_PC[200], apPC[200]
for i = 0,199 {
PCcell[i] apc_PC[i] = new APCount(0.5)
apc_PC[i].thresh = -45
apPC[i] = new Vector()
apc_PC[i].record(apPC[i])
}

// Record ION
objref apc_ION[40], apION[40]
for i = 0,39 {
IONcell[i] apc_ION[i] = new APCount(0.5)
apc_ION[i].thresh = -45
apION[i] = new Vector()
apc_ION[i].record(apION[i])
}

//////////////////////////////////////////////////////////////////
// Simulation starts											//
//																//
// Set up synaptic transmission delays between PC-DCN and PC-NO	//
continuerun(t+tstop)											//
//																//
//////////////////////////////////////////////////////////////////

// Record DCN
objref rec_v_DCN1, rec_v_DCN2, rec_v_DCN3, rec_v_DCN4, rec_v_DCN5
rec_v_DCN1 = new Vector()
rec_v_DCN1.append(DCNcell[0].v(0.5))
rec_v_DCN2 = new Vector()
rec_v_DCN2.append(DCNcell[1].v(0.5))
rec_v_DCN3 = new Vector()
rec_v_DCN3.append(DCNcell[2].v(0.5))
rec_v_DCN4 = new Vector()
rec_v_DCN4.append(DCNcell[3].v(0.5))
rec_v_DCN5 = new Vector()
rec_v_DCN5.append(DCNcell[4].v(0.5))

// Record PC
objref rec_v_PC
rec_v_PC = new Vector()
rec_v_PC.append(PCcell[140].v(0.5))

// Record Vim
objref rec_v_Vim1, rec_v_Vim2, rec_v_Vim3, rec_v_Vim4, rec_v_Vim5
rec_v_Vim1 = new Vector()
rec_v_Vim1.append(TCcell[0].v(0.5))
rec_v_Vim2 = new Vector()
rec_v_Vim2.append(TCcell[1].v(0.5))
rec_v_Vim3 = new Vector()
rec_v_Vim3.append(TCcell[2].v(0.5))
rec_v_Vim4 = new Vector()
rec_v_Vim4.append(TCcell[3].v(0.5))
rec_v_Vim5 = new Vector()
rec_v_Vim5.append(TCcell[4].v(0.5))

// Record SIN
objref rec_sinuStim
rec_sinuStim = new Vector()
rec_sinuStim.append(sinuStim_group[0].i)

// Save PC
objref rec_PC
rec_PC = new File()
rec_PC.aopen(f_PCap)
for i = 0,199 {
    if (apPC[i].size()>0) {rec_PC.printf("%d,%f\n",i,apPC[i].x[0])}
}
rec_PC.close()

// Save DCN
objref sav_v_DCN1, sav_v_DCN2, sav_v_DCN3, sav_v_DCN4, sav_v_DCN5
sav_v_DCN1 = new File()
sav_v_DCN1.aopen(f_DCN1)
rec_v_DCN1.printf(sav_v_DCN1)
sav_v_DCN1.close()

sav_v_DCN2 = new File()
sav_v_DCN2.aopen(f_DCN2)
rec_v_DCN2.printf(sav_v_DCN2)
sav_v_DCN2.close()

sav_v_DCN3 = new File()
sav_v_DCN3.aopen(f_DCN3)
rec_v_DCN3.printf(sav_v_DCN3)
sav_v_DCN3.close()

sav_v_DCN4 = new File()
sav_v_DCN4.aopen(f_DCN4)
rec_v_DCN4.printf(sav_v_DCN4)
sav_v_DCN4.close()

sav_v_DCN5 = new File()
sav_v_DCN5.aopen(f_DCN5)
rec_v_DCN5.printf(sav_v_DCN5)
sav_v_DCN5.close()

// Save PC
objref sav_v_PC
sav_v_PC = new File()
sav_v_PC.aopen(f_PCv)
rec_v_PC.printf(sav_v_PC)
sav_v_PC.close()

// Save Vim
objref sav_v_Vim1, sav_v_Vim2, sav_v_Vim3, sav_v_Vim4, sav_v_Vim5
sav_v_Vim1 = new File()
sav_v_Vim1.wopen(f_Vim1)
rec_v_Vim1.printf(sav_v_Vim1)
sav_v_Vim1.close()

sav_v_Vim2 = new File()
sav_v_Vim2.wopen(f_Vim2)
rec_v_Vim2.printf(sav_v_Vim2)
sav_v_Vim2.close()

sav_v_Vim3 = new File()
sav_v_Vim3.wopen(f_Vim3)
rec_v_Vim3.printf(sav_v_Vim3)
sav_v_Vim3.close()

sav_v_Vim4 = new File()
sav_v_Vim4.wopen(f_Vim4)
rec_v_Vim4.printf(sav_v_Vim4)
sav_v_Vim4.close()

sav_v_Vim5 = new File()
sav_v_Vim5.wopen(f_Vim5)
rec_v_Vim5.printf(sav_v_Vim5)
sav_v_Vim5.close()

// Save ION
objref rec_ION
rec_ION = new File()
rec_ION.aopen(f_ION)
for i = 0,39 {
    if (apION[i].size()>0) {rec_ION.printf("%d,%f\n",i,apION[i].x[0])}
}
rec_ION.close()

// Save SIN
objref sav_sinuStim
sav_sinuStim = new File()
sav_sinuStim.aopen(f_sin)
rec_sinuStim.printf(sav_sinuStim)
sav_sinuStim.close()

objref ss,f_ss
ss = new SaveState()
f_ss = new File()
ss.save()
f_ss.wopen("currentState.dat")
ss.fwrite(f_ss)
f_ss.close()

quit()