////////////////////
/* net_dd_ana.hoc */
////////////////////
/* -------------------------------------------------------------------------- */
/* -------------- ANALYSIS -------------------------------------------------- */
/* -------------------------------------------------------------------------- */
strdef strvar // str variable for file name template
objref dfile
dfile = new File()
objref AP_IN_cells, AP_IN_times, AP_PN_cells, AP_PN_times, Cellvar_IN, Cellvar_PN
objref AP_IN_mx[nIN+1], AP_PN_mx[nPN+1] // matrix for binned AP data
Cellvar_IN = new Vector(nIN)
Cellvar_PN = new Vector(nPN)
binw = 0.2 // bin width for evaluating coherence - up to 500 Hz
bins = int(tstop/binw) // n of bins
/* Saving of the APs ---------------------------------------------------------*/
AP_IN_cells = new Vector((tstop*500/1000)*nIN) // two vectors to store cell_i & time of APs
AP_IN_times = new Vector((tstop*500/1000)*nIN) // assumes a max firing freq of 500Hz!
AP_IN_count = 0
AP_PN_cells = new Vector((tstop*500/1000)*nPN)
AP_PN_times = new Vector((tstop*500/1000)*nPN)
AP_PN_count = 0
for i = 0,nIN {
AP_IN_mx[i] = new Vector(bins+1) // array of vectors. For each cell, there is one
} // vector which has for each time-step-bin one element
// in which a "1" indicates an AP and a "0" indicates no AP.
for i = 0,nPN {
AP_PN_mx[i] = new Vector(bins+1)
}
proc setbin() {local b
b = $1
if (b<0.2) {b = 0.2} // do not allow smaller values then 0.2 ms for bin size
// this would mean an >500Hz osc. we don't look at that now
binw = b // matrix was sized to this max.
bins = int(tstop/binw)
}
proc AP2_mx() {local i
// form a bin matrix from AP time data vectors
for i = 0,nIN { AP_IN_mx[i].fill(0) }
for i = 0,AP_IN_count-1 {
AP_IN_mx[AP_IN_cells.x[i]].x[int(AP_IN_times.x[i]/binw)] = 1
AP_IN_mx[nIN].x[int(AP_IN_times.x[i]/binw)] += 1
// count number of active cells
// for the histogram
}
for i = 0,nPN { AP_PN_mx[i].fill(0) }
for i = 0,AP_PN_count-1 {
AP_PN_mx[AP_PN_cells.x[i]].x[int(AP_PN_times.x[i]/binw)] = 1
AP_PN_mx[nPN].x[int(AP_PN_times.x[i]/binw)] += 1
// count number of active cells
// for the histogram
}
}
proc saveHist() { local i
// interneurons
sprint(strvar,"%si_hist.tab",filestr)
dfile.wopen(strvar)
setbin(1)
AP2_mx() // bin data at 1 ms
for i = int(hbcoh_start)/binw,int(hbcoh_end)/binw {
dfile.printf("%f\t%d\n",i*binw,AP_IN_mx[nIN].x[i]) // no of active cells stored here
}
dfile.close()
rtime = stopsw()
printf("%d:%d -- icell-Histogram data saved in file %s. \n",rtime/60,rtime%60,strvar)
// principal neurons
sprint(strvar,"%se_hist.tab",filestr)
dfile.wopen(strvar)
for j = int(hbcoh_start)/binw,int(hbcoh_end)/binw {
dfile.printf("%f\t%d\n",j*binw,AP_PN_mx[nPN].x[j]) // no of active cells stored here
}
dfile.close()
rtime = stopsw()
printf("%d:%d -- ecell-Histogram data saved in file %s. \n",rtime/60,rtime%60,strvar)
}
proc save_raster() { // Saves the spike data in order to rebuild the raster plot.
// interneurons
sprint(strvar,"%si_raster.tab",filestr)
dfile.wopen(strvar)
dfile.printf("Time\tCell no\n")
for i = 0,AP_IN_count-1 { dfile.printf("%f\t%d\n",AP_IN_times.x[i],AP_IN_cells.x[i]) }
dfile.close()
// principal neurons
sprint(strvar,"%se_raster.tab",filestr)
dfile.wopen(strvar)
dfile.printf("Time\tCell no\n")
for i = 0,AP_PN_count-1 { dfile.printf("%f\t%d\n",AP_PN_times.x[i],AP_PN_cells.x[i]) }
dfile.close()
}
/* Saving of the synaptic conductances ---------------------------------------*/
objref gII, gIE, gEI
if (export_traces == 1) {
gII = new Matrix(nIN,1+tstop/t_step)
gIE = new Matrix(nPN,1+tstop/t_step)
gEI = new Matrix(nIN,1+tstop/t_step)
}
func get_total_gII() { local s,j
s = 0
for j = 0,Syn_II_N-1 {
s = s + interneuron[$1].syn_II[j].g
}
return s
}
proc write_gII() { local i // $1 time step --> element of the record matrix
for i = 0,nIN-1 {
gII.x[i][$1] = get_total_gII(i)
}
}
proc save_gII() {
sprint(strvar,"%sgII.tab",filestr)
dfile.wopen(strvar)
gII.fprint(0,dfile)
}
func get_total_gIE() { local s,j
s = 0
for j = 0,Syn_IE_N-1 {
s = s + principalneuron[$1].syn_IE[j].g
}
return s
}
proc write_gIE() { local i // $1 time step --> element of the record matrix
for i = 0,nPN-1 {
gIE.x[i][$1] = get_total_gIE(i)
}
}
proc save_gIE() {
sprint(strvar,"%sgIE.tab",filestr)
dfile.wopen(strvar)
gIE.fprint(0,dfile)
}
proc write_gEI() { local i // $1 time step --> element of the record matrix
for i = 0,nIN-1 {
gEI.x[i][$1] = interneuron[i].syn_EI.g
}
}
proc save_gEI() {
sprint(strvar,"%sgEI.tab",filestr)
dfile.wopen(strvar)
gEI.fprint(0,dfile)
}
/* Saving of the inhibitory synaptic currents (LFP) --------------------------*/
objref lfp
if (record_lfp == 1) {
lfp = new Matrix(nIN,1+tstop/t_step)
}
func get_moment_lfp() { local j,s,pirat,p // $1: position within network (#IN)
s = 0
for j = 0,Syn_II_N-1 {
s = s + interneuron[$1].syn_II[j].g
}
for pirat = 0,int(nPN/nIN)-1 {
p = ($1*int(nPN/nIN))+pirat
for j = 0,Syn_IE_N-1 {
s = s + principalneuron[p].syn_IE[j].g
}
}
return s
}
proc save_lfp() {
sprint(strvar,"%slfp.tab",filestr)
dfile.wopen(strvar)
lfp.fprint(0,dfile)
}
/* Saving the connectivity matrices ----------------------------------------- */
objref conmx_II, conmx_IE, conmx_EI, conmx_EE
proc save_conmx() {
sprint(strvar,"%sconmxII.tab",filestr)
dfile.wopen(strvar)
conmx_II.fprint(dfile)
sprint(strvar,"%sconmxIE.tab",filestr)
dfile.wopen(strvar)
conmx_IE.fprint(dfile)
sprint(strvar,"%sconmxEI.tab",filestr)
dfile.wopen(strvar)
conmx_EI.fprint(dfile)
sprint(strvar,"%sconmxEE.tab",filestr)
dfile.wopen(strvar)
conmx_EE.fprint(dfile)
}
func ind_gcomp() { local j, G // $1 gsyn, $2 rise, $3 decay
for j = 0,(100/dt)-1 {
G = G+($1*(exp(-(j*dt)/$3)-exp(-(j*dt)/$2)))
}
return G
}
/* Saving of the membrane voltages -------------------------------------------*/
objref vIN, vPN
if (export_traces == 1) {
vIN = new Matrix(nIN,1+tstop/t_step)
vPN = new Matrix(nPN,1+tstop/t_step)
}
proc write_vIN() { local i // $1 time step --> element of the record matrix
for i = 0,nIN-1 {
vIN.x[i][$1] = interneuron[i].soma_interneuron.v(0.5)
}
}
proc write_vPN() { local i // $1 time step --> element of the record matrix
for i = 0,nPN-1 {
vPN.x[i][$1] = principalneuron[i].soma_principalneuron.v(0.5)
}
}
proc save_vIN() {
sprint(strvar,"%svIN.tab",filestr)
dfile.wopen(strvar)
vIN.fprint(0,dfile)
}
proc save_vPN() {
sprint(strvar,"%svPN.tab",filestr)
dfile.wopen(strvar)
vPN.fprint(0,dfile)
}
/* -------------------------------------------------------------------------- */
/* -------------- PLOTTING -------------------------------------------------- */
/* -------------------------------------------------------------------------- */
objref v_IN_plt, v_PN_plt
proc initv_IN_Plt() { local i // plot voltage traces
v_IN_plt = new Graph(0) // init Graph object instance
v_IN_plt.size(tsyn,tstop,-80,80) // window size
v_IN_plt.label(0.95, 0.58, "ms") // x label
v_IN_plt.label(0.01, 0.95, "nS") // y label
v_IN_plt.label(0.01, 0.80, "EPSC") // y label
v_IN_plt.label(0.01, 0.65, "IPSC") // y label
v_IN_plt.label(0.01, 0.55, "mV") // y label
v_IN_plt.label(0.95, 0.9, "IN")
v_IN_plt.label(1,0) // move for next labels out
v_IN_plt.view(tsyn, -80, tstop-tsyn, 160, 0, 20, 300, 230)
for i = 0,mx_v_IN-1 {
sprint(strvar,"%s%d%s","interneuron[n_v_IN_",i,"].soma_interneuron.v(0.5)")
v_IN_plt.addvar(strvar,i%4+1 ,1)
}
}
proc initv_PN_Plt() { local i // plot voltage traces
v_PN_plt = new Graph(0) // init Graph object instance
v_PN_plt.size(tsyn,tstop,-80,80) // window size
v_PN_plt.label(0.95, 0.58, "ms") // x label
v_PN_plt.label(0.01, 0.95, "nS") // y label
v_PN_plt.label(0.01, 0.80, "EPSC") // y label
v_PN_plt.label(0.01, 0.65, "IPSC") // y label
v_PN_plt.label(0.01, 0.55, "mV") // y label
v_PN_plt.label(0.95, 0.9, "PN")
v_PN_plt.label(1,0) // move for next labels out
v_PN_plt.view(tsyn, -80, tstop-tsyn, 160, 0, 320, 300, 230)
for i = 0,mx_v_PN-1 {
sprint(strvar,"%s%d%s","principalneuron[n_v_PN_",i,"].soma_principalneuron.v(0.5)")
v_PN_plt.addvar(strvar,i%4+1 ,1)
}
}
objref r_plt, h_plt
proc initr_Plt() {
r_plt = new Graph(0) // init Graph object instance
r_plt.size(0,tstop,0,nPN) // window size
r_plt.label(0.95, 0.02, "ms") // x label
r_plt.label(0.01, 0.81, "cell#") // y label
r_plt.color(2)
r_plt.label(0.9, 0.9, "INs") // cell type
r_plt.color(4)
r_plt.label("PNs") // cell type
r_plt.view(0,0,tstop, nPN, 320, 20, 300, 230)
}
proc inithPlt() { // plot raster of APs (cell x time matrix)
h_plt = new Graph(0) // init Graph object instance
h_plt.size(0,tstop,0,0.5) // window size
h_plt.label(0.95, 0.02, "ms") // x label
h_plt.label(0.01, 0.81, "% cell") // y label
h_plt.color(2)
h_plt.label(0.9, 0.9, "INs") // cell type
h_plt.color(4)
h_plt.label("PNs") // cell type
h_plt.view(0,0,tstop,0.5, 640, 20, 300, 230)
}
proc plotAPs() {
r_plt.erase()
for i = 0,AP_IN_count-1 {
r_plt.mark(AP_IN_times.x[i],(AP_IN_cells.x[i]+1)*(nPN/nIN),"o",1,2,1)
}
for j = 0,AP_PN_count-1 {
r_plt.mark(AP_PN_times.x[j],AP_PN_cells.x[j]+1,"o",1,4,1)
}
r_plt.flush()
doEvents()
}
proc plotHist() { local i
setbin(1) AP2_mx() // we need the matrix first
h_plt.erase()
h_plt.beginline(2,1)
for i = 0,bins-1 { h_plt.line(i*binw,AP_IN_mx[nIN].x[i]/nIN) }
h_plt.beginline(4,1)
for i = 0,bins-1 { h_plt.line(i*binw,AP_PN_mx[nPN].x[i]/nPN) }
h_plt.flush()
doEvents()
}
proc saveAPPlt(){
sprint(strvar,"%sAP_plot.eps",filestr)
r_plt.printfile(strvar)
rtime = stopsw()
printf("%d:%d -- Raster plot saved in file %s. \n",rtime/60,rtime%60,strvar)
}
proc savehPlt(){
sprint(strvar,"%sh_plot.eps",filestr)
h_plt.printfile(strvar)
}
