initialization_tstart = startsw()
strdef modelname
strdef cmd
strdef fileName
objref fileObj
objref fileObj2
fileObj = new File()
fileObj2 = new File()
rdSeed = 1
Nmc = 120 //the number of PCs should be 110,120,130 in different conditions, corresponding to Fig 6C-E
trial =1
con_rad = 5
beta = 0
isolat =0
pre_pf_rec = 0
pre_inh_rec =0
celsius = 34
tstop = 500000
dt = 0.02
step_per_ms= 1/dt
v_init = -70
NsynI = 10
Econ = 0.0005
Icon = 0.001/NsynI
NcontE = 2
NcontI = 5
Ncontst = 6*NsynI
Ncontbs = 4*NsynI
frpf_scale = 1e-6
frin_scale = 1e-6
Ihold = 0
Ivar = 0.1940
param_NGrc_sign = 0
param_Nst_sign = 0
param_Nbs_sign = 0
NGrc = 1
Nst = 1
Nbs = 1
param_f_s = 2
param_f_s2 = 1
param_amp_s = 0
param_f_inh_s = 2
param_f_inh_s2 = 1
param_amp_inh_s = 0
strdef outDir,outDir2,outDir3
sprint(outDir,"predata")
sprint(cmd,"system(\"mkdir -p %s\")",outDir)
execute(cmd)
NsynE = (NGrc+param_NGrc_sign)*NcontE
NsynIst = (Nst+param_Nst_sign)*Ncontst
NsynIbs = (Nbs+param_Nbs_sign)*Ncontbs
param_rateE =0.5400*frpf_scale
param_rateIst = 14.4*frin_scale
param_rateIbs = 14.4*frin_scale
xopen("pkj.hoc")
objref tempvec
tempvec = new Vector()
{tempvec.append(Nmc)}
{Ncells2save = tempvec.min()}
objref pkj_TTC
objref sl //synaptic locations list
objref rds1,rds2
{rds1 = new Random(1000*rdSeed)}
{rds2 = new Random(2000*rdSeed)}
{rds1.uniform(0,1)}
{rds2.uniform(0,1)}
objref conMat,cons,cont,vecN,switchEdge,target,conInd,switchInd,switchInd2,si
conMat = new Matrix(Nmc,Nmc)
cons = new Matrix(Nmc,con_rad)
cont = new Matrix(Nmc,con_rad)
vecN = new Vector()
vecN.indgen(0,Nmc-1,1)
for (i=0;i<con_rad;i+=1) {
cons.setcol(i,vecN)
cont.setcol(i,vecN.c.add(i+1))
}
for(i=0;i<Nmc;i+=1){
for(j=0;j<con_rad;j+=1){
cont.x[i][j] = (cont.x[i][j])%Nmc
conMat.x[cons.x[i][1]][cont.x[i][j]] = 1
}
}
conMat.add(conMat.c.transpose())
switchEdge = new Vector(2*con_rad)
switchInd = new Vector()
switchInd2= new Vector()
target = new Vector(Nmc)
conInd = new Vector()
for (i = 0;i<Nmc;i+=1) {
switchEdge.setrand(rds1)
target.setrand(rds2)
target.x(i) = 0
switchInd = switchEdge.c.indvwhere("<",beta)
switchInd2 = switchEdge.c.indvwhere(">=",beta)
conInd = conMat.getcol(i).c.indvwhere(">",0)
if (switchInd.size()*conInd.size()) {
for (j = 0;j<switchInd2.size();j+=1) {
target.x(conInd.x(switchInd2.x(j))) = 0
}
si = target.sortindex.reverse()
for (k = 0;k<switchInd.size();k+=1) {
conMat.x[conInd.x(switchInd.x[k])][i]=0
}
for (k = 0;k<switchInd.size();k+=1) {
conMat.x[si.x[k]][i]=1
}
}
}
objref het_vec,hetrand,het_index
hetrand = new Random(6666666+Nmc)
hetrand.uniform(0,1)
het_vec = new Vector(Nmc)
het_vec.setrand(hetrand)
het_index = het_vec.sortindex
if (isolat) {
conMat.zero()
con_rad = 0
}
{load_file("netparmpi.hoc")}
objref epnm
epnm = new ParallelNetManager(Nmc)
{epnm.round_robin()}
for(i=0;i<Nmc;i+=1){
if (epnm.gid_exists(i)) {
pkj_TTC = new pkj_neuron()
epnm.register_cell(i,pkj_TTC)
epnm.pc.gid2cell(i).initRand(1000*rdSeed+i)
epnm.pc.gid2cell(i).setnetworkparameters(Econ,Icon,NsynE,NsynIst,NsynIbs,NcontI)
epnm.pc.gid2cell(i).setsignalparameters(param_rateE,param_rateIst,NGrc,Nst, Nbs, param_NGrc_sign,param_Nst_sign,param_Nbs_sign,param_f_s,param_amp_s,param_f_inh_s,param_amp_inh_s,tstop)
epnm.pc.gid2cell(i).distributeSyn()
epnm.pc.gid2cell(i).setpretrains(i,i+500,i+1000,i+1500,i+2000)
}
}
proc distribute() {
Nred = Nmc-100
het_index.resize(Nred)
for(i=0;i<Nmc;i+=1){
if (epnm.gid_exists(i)) {
epnm.pc.gid2cell(i).fl.new_seed(i+trial*Nmc)
epnm.pc.gid2cell(i).fl.m = -0.2+($1-1)*0.1
epnm.pc.gid2cell(i).fl.s = 0+Ivar*($2-1)*0.2
epnm.pc.gid2cell(i).fl.tau = 5
epnm.pc.gid2cell(i).fl.amp = 1
if (het_index.contains(i)) {
epnm.pc.gid2cell(i).fl.amp = -0.2
}
epnm.pc.gid2cell(i).fl.freq = 1
epnm.pc.gid2cell(i).fl.fr_stm = 1
epnm.pc.gid2cell(i).fl.onset = 0
epnm.pc.gid2cell(i).fl.np = 100000000
epnm.pc.gid2cell(i).Istim.dur = i*3
}
}
}
objref syninds
proc connection() {
for(i=0;i<Nmc;i+=1){
if (epnm.gid_exists(i)) {
epnm.pc.gid2cell(i).insertMCcons(conMat.getcol(i))
// distribute synapses on the postsynaptic cell.
}
}
syninds = new Vector()
for(i=0;i<Nmc;i+=1){
syninds.append(NsynE+NsynIst+NsynIbs) //syninds starts from the total number of background synapses and then ....
}
for(i=0;i<Nmc;i+=1){
for(j=0;j<Nmc;j+=1){
if (conMat.x[j][i] != 0){
for(jj=0;jj<NcontI;jj+=1){
epnm.nc_append(j,i,syninds.x[i],1e-3/NcontI,1.5)
syninds.x[i] +=1
}
}
}
}
for(i=0;i<Nmc;i+=1){
if (epnm.gid_exists(i)) {
{epnm.pc.gid2cell(i).queuePreTrains()}
}1
}
}
connection()
strdef outDir2
proc out_folder() {
modelname = "Purkinje"
sprint(outDir2,"simdata_sango/%s/N%d_rad%d_beta%d_stimm%d_stims%d_ratio%d",modelname,Nmc,con_rad,beta*10,$1,$2,$3)
sprint(cmd,"system(\"mkdir -p %s\")",outDir2)
execute(cmd)
if(epnm.gid_exists(0)){
sprint(fileName,"%s/conmat.txt",outDir2)
fileObj.wopen(fileName)
conMat.fprint(0,fileObj)
fileObj.close()
sprint(fileName,"%s/het_rand.txt",outDir2)
fileObj.wopen(fileName)
het_vec.printf(fileObj)
fileObj.close()
}
}
objref vSomaList
objref apcvecList, apcList
proc record() {
apcvecList = new List()
apcList = new List()
vSomaList = new List()
for(i=0;i<Ncells2save;i+=1){
if (epnm.gid_exists(i)) {
vSomaList.append(new Vector())
access epnm.pc.gid2cell(i).somaA
vSomaList.o[vSomaList.count()-1].record(&v(0.5),dt)
}
}
for(i=0;i<Nmc;i+=1){
if (epnm.gid_exists(i)) {
access epnm.pc.gid2cell(i).somaA
apcList.append(new APCount(0.5))
apcvecList.append(new Vector())
apcList.o[apcList.count()-1].thresh= -30
apcList.o[apcList.count()-1].record(apcvecList.o[apcList.count()-1])
}
}
}
proc myprun() {
{epnm.set_maxstep(100)}
finitialize(-65)
if (epnm.gid_exists(0)) {
print "\n"
sim_tstart = startsw()
initializationtime = (sim_tstart-initialization_tstart)/3600
print "Initialization completed. Initialization took ", initializationtime, " hours\n"
print "Starting simulation\n"
print "\n"
}
{epnm.psolve(tstop)}
if (epnm.gid_exists(0)) {
simruntime = (startsw() - sim_tstart)/3600
print "Simulation took ", simruntime, " hours\n"
fileObj2.printf("simulation took %f hours\n",simruntime)
fileObj2.flush()
}
if (pre_inh_rec) {
for(i=0;i<Ncells2save;i+=1){
if (epnm.gid_exists(i)) {
for (mm=0; mm< epnm.pc.gid2cell(i).preTrainList_inh_sign.count();mm+=1) {
tempvec = new Vector()
tempvec = epnm.pc.gid2cell(i).preTrainList_inh_sign.o[mm].c
sprint(fileName,"%s/pre_inh_%03d_%03d.txt",outDir,i,mm)
fileObj.wopen(fileName)
tempvec.printf(fileObj,"%2.6e\n")
fileObj.close()
}
}
}
}
if (pre_pf_rec) {
for(i=0;i<Ncells2save;i+=1){
if (epnm.gid_exists(i)) {
for (mm=0; mm< epnm.pc.gid2cell(i).preTrainList_sign.count();mm+=1) {
tempvec = new Vector()
tempvec = epnm.pc.gid2cell(i).preTrainList_sign.o[mm].c
sprint(fileName,"%s/pre_pf_%03d_%03d.txt",outDir,i,mm)
fileObj.wopen(fileName)
tempvec.printf(fileObj,"%2.6e\n")
fileObj.close()
}
}
}
}
i2=0
for(i=0;i<Nmc;i+=1){
if (epnm.gid_exists(i)) {
tempvec = new Vector()
tempvec = apcvecList.o[i2].c
sprint(fileName,"%s/spikeTimesN%d.txt",outDir2,i)
fileObj.wopen(fileName)
tempvec.printf(fileObj,"%2.6e\n")
fileObj.close()
tempvec = new Vector()
Nspike = apcList.o[i2].n
i2+=1
}
}
}
proc loop_run() {
for iloop =3,3 {
for jloop = 7,7 {
for mloop = 2,2 {
distribute(iloop,jloop,mloop)
out_folder(iloop,jloop,mloop)
record()
myprun()
}
}
}
}
sprint(outDir3,"simtime")
sprint(cmd,"system(\"mkdir -p %s\")",outDir3)
execute(cmd)
fileObj2.wopen("./simtime/simulationtime.dat")
loop_run()
fileObj2.close()
{epnm.pc.runworker()}
{epnm.pc.done()}