// calcisilag.hoc // Olfactory bulb network model: calculate inter-spike interval // : and lag time statistics // Andrew Davison, The Babraham Institute, 2000. /* The following procedures are defined for writing results to file * fileroot is the filename root. A suffix will be added to * this, e.g. fileroot.synch for print_si(). * i,j are the mitral cell indices * * print_smooth_hist(variance,fileroot) * print_gran_smooth_hist(variance,fileroot) * print_spiketimes(i,j,fileroot) * print_raster(fileroot) * print_gran_raster(fileroot) * print_isis(i,j,fileroot) * print_isi_stats(fileroot) * print_lags(i,j,fileroot) * print_si(fileroot) */ // Variables used in this file objref work, work2, outputarray, isi, lags, hist work = new Vector() work2 = new Vector() outputarray = new Matrix(nmitx,nmity) isi = new Vector() lags = new Vector() // Procedures for processing spike times ------------------------------- proc calc_isis() { local i,j,k,n // 3 args - indices of mitral cell, transient time if (\$1 > nmitx || \$2 > nmity) { print "Sorry - index out of range. Please try again." return } i = int(\$1) j = int(\$2) isi.resize(0) n = mit[i][j].spiketimes.size() if (n > 1) { for k = 1,n-1 { if (mit[i][j].spiketimes.x[k-1] > \$3) { isi.append(mit[i][j].spiketimes.x[k]-mit[i][j].spiketimes.x[k-1]) } } } } proc calc_gran_isis() { local i,j,k,n // 3 args - indices of granule cell, transient time if (\$1 > ngranx || \$2 > ngrany) { print "Sorry - index out of range. Please try again." return } i = int(\$1) j = int(\$2) isi.resize(0) n = gran[i][j].spiketimes.size() if (n > 1) { for k = 1,n-1 { if (gran[i][j].spiketimes.x[k-1] > \$3) { isi.append(gran[i][j].spiketimes.x[k]-gran[i][j].spiketimes.x[k-1]) } } } } func minisi() { local i,j,min // find shortest mean ISI min = 1e6 for i = 0, nmitx-1 { for j = 0, nmity-1 { calc_isis(i,j,ttrans) if (isi.size() > 0) { if (isi.mean() < min) { min = isi.mean() } } } } return min } proc rate_array() { local i,j for i = 0,nmitx-1 { for j = 0, nmity-1 { calc_isis(i,j,ttrans) if (isi.size() > 0) { outputarray.x[i][j] = 1000/isi.mean() } else { outputarray.x[i][j] = 0 } } } print "Max: ",arraymax(outputarray) outputarray.muls(1/arraymax(outputarray)) } proc calc_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of mitral cells, transient time if (\$1 > nmitx || \$2 > nmity || \$3 > nmitx || \$4 > nmity) { print "Sorry - index out of range. Please try again." return } i1 = int(\$1) j1 = int(\$2) i2 = int(\$3) j2 = int(\$4) lags.resize(0) // for each spiketime in cell 1, find closest spike in cell 2 // Note: first and last spikes ignored since can't calculate previous ISI if (mit[i2][j2].spiketimes.size > 0) { for k = 1,mit[i1][j1].spiketimes.size()-2 { if (mit[i1][j1].spiketimes.x[k] > \$5) { work = mit[i2][j2].spiketimes.c.add(-mit[i1][j1].spiketimes.x[k]) minidx = work.c.abs.min_ind() min = work.x[minidx] isiprev = mit[i1][j1].spiketimes.x[k-1]-mit[i1][j1].spiketimes.x[k] isinext = mit[i1][j1].spiketimes.x[k+1]-mit[i1][j1].spiketimes.x[k] if (min > isiprev/2 && min < isinext/2) { lags.append(min) } } } } } proc calc_phase_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of mitral cells, transient time if (\$1 > nmitx || \$2 > nmity || \$3 > nmitx || \$4 > nmity) { print "Sorry - index out of range. Please try again." return } i1 = int(\$1) j1 = int(\$2) i2 = int(\$3) j2 = int(\$4) lags.resize(0) // for each spiketime in cell 1, find closest spike in cell 2 // Note: first and last spikes ignored since can't calculate previous ISI if (mit[i2][j2].spiketimes.size > 0) { for k = 1,mit[i1][j1].spiketimes.size()-2 { if (mit[i1][j1].spiketimes.x[k] > \$5) { work = mit[i2][j2].spiketimes.c.add(-mit[i1][j1].spiketimes.x[k]) minidx = work.c.abs.min_ind() min = work.x[minidx] isiprev = mit[i1][j1].spiketimes.x[k-1]-mit[i1][j1].spiketimes.x[k] isinext = mit[i1][j1].spiketimes.x[k+1]-mit[i1][j1].spiketimes.x[k] if (min > isiprev/2 && min < isinext/2) { if (min < 0) { lags.append(min/isiprev) } else { lags.append(min/isinext) } } } } } } proc calc_gran_lags() { local i1,j1,i2,j2,k,minidx,min // 5 args - indices of granule cells, transient time if (\$1 > ngranx || \$2 > ngrany || \$3 > ngranx || \$4 > ngrany) { print "Sorry - index out of range. Please try again." return } i1 = int(\$1) j1 = int(\$2) i2 = int(\$3) j2 = int(\$4) lags.resize(0) // for each spiketime in cell 1, find closest spike in cell 2 // Note: first and last spikes ignored since can't calculate previous ISI if (gran[i2][j2].spiketimes.size > 0) { for k = 1,gran[i1][j1].spiketimes.size()-2 { if (gran[i1][j1].spiketimes.x[k] > \$5) { work = gran[i2][j2].spiketimes.c.add(-gran[i1][j1].spiketimes.x[k]) minidx = work.c.abs.min_ind() min = work.x[minidx] isiprev = gran[i1][j1].spiketimes.x[k]-gran[i1][j1].spiketimes.x[k-1] isinext = gran[i1][j1].spiketimes.x[k]-gran[i1][j1].spiketimes.x[k+1] if (min < isiprev/2 && min > isinext/2) { lags.append(min) } } } } } proc time_hist() { // 1 arg - time step work.resize(0) for i = 0, nmitx-1 { for j = 0, nmity-1 { work.append(mit[i][j].spiketimes) } } hist = work.histogram(0,tstop,\$1) hist.printf("%d\n") } func synch_index() { local i1,j1,i2,j2,n synchindex = 0 n = 0 for i1 = 0, nmitx-1 { for j1 = 0, nmity-1 { if (mit[i1][j1].spiketimes.size() > 0) { for i2 = 0, nmitx-1 { for j2 = 0, nmity-1 { if (i1 != i2 || j1 != j2) { calc_phase_lags(i1,j1,i2,j2,ttrans) n += lags.size() synchindex += lags.reduce("abs",0) } } } } } } if (n > 0) { return synchindex/n } else { return 1e6 } } func phaselock_index() { local n,i1,j1,i2,j2 synchindex = 0 n = 0 for i1 = 0, nmitx-1 { for j1 = 0, nmity-1 { if (mit[i1][j1].spiketimes.size() > 0) { for i2 = 0, nmitx-1 { for j2 = 0, nmity-1 { if (i1 != i2 || j1 != j2) { calc_phase_lags(i1,j1,i2,j2,ttrans) if (lags.size() > 1) { synchindex += lags.var() n += 1 } } } } } } } if (n > 0) { synchindex = sqrt(synchindex/n) return synchindex } else { return 1e6 } } // Procedures for writing out data -------------------------------- proc calc_smooth_hist() { // 1 arg - variance work.resize(0) for i = 0, nmitx-1 { for j = 0, nmity-1 { work.append(mit[i][j].spiketimes) } } hist = work.sumgauss(0,tstop,1,\$1) } proc print_smooth_hist() { // 2 args - variance, filename root calc_smooth_hist(\$1) sprint(filename,"%s.smhist",\$s2) outfile.wopen(filename) outfile.printf("# Mitral cell smoothed histogram\n") hist.printf(outfile,"%8.3f\n") outfile.close() /* work.resize(0) hist.remove(0,ttrans) work.spctrm(hist) sprint(filename,"%s.pow",\$s2) outfile.wopen(filename) outfile.printf("# Power spectrum of Mitral cell smoothed histogram\n") work.printf(outfile,"%9.5f\n") outfile.close() */ } proc calc_gran_smooth_hist() { // 1 arg - variance work.resize(0) for i = 0, ngranx-1 { for j = 0, ngrany-1 { work.append(gran[i][j].spiketimes) } } hist = work.sumgauss(0,tstop,1,\$1) } proc print_gran_smooth_hist() { // 2 args - variance, filename root calc_gran_smooth_hist(\$1) sprint(filename,"%s.gran.smhist",\$s2) outfile.wopen(filename) outfile.printf("# Granule cell smoothed histogram\n") hist.printf(outfile,"%8.3f\n") outfile.close() //work.resize(0) //hist.remove(0,ttrans) //work.spctrm(hist) //sprint(filename,"%s.gran.pow",\$s2) //outfile.wopen(filename) //outfile.printf("# Power spectrum of Granule cell smoothed histogram\n") //work.printf(outfile,"%9.5f\n") //outfile.close() } proc print_gran_hist() { // 2 args - binsize, filename root work.resize(0) for i = 0, ngranx-1 { for j = 0, ngrany-1 { work.append(gran[i][j].spiketimes) } } hist = work.histogram(0,tstop,\$1) sprint(filename,"%s.gran.hist",\$s2) outfile.wopen(filename) outfile.printf("# Granule cell unsmoothed histogram\n") hist.printf(outfile,"%8.3f\n") outfile.close() } proc print_spiketimes() { // 3 args - indices of mitral cell plus filename root if (numarg() == 3) { sprint(filename,"%s_%d_%d.isi",\$s3,\$1,\$2) outfile.wopen(filename) outfile.printf("# Spiketimes for mitral cell [%d][%d]",\$1,\$2) mit[\$1][\$2].spiketimes.printf(outfile,"%10.3f") } if (numarg() == 1) { sprint(filename,"%s.isi",\$s1) outfile.wopen(filename) for i = 0, nmitx-1 { for j = 0, nmity-1 { outfile.printf("[%d][%d]",i,j) mit[i][j].spiketimes.printf(outfile,"%10.3f") } } } outfile.close() } proc print_raster() { // 1 arg - filename root sprint(filename,"%s.ras",\$s1) outfile.wopen(filename) outfile.printf("# Mitral cell raster plot\n") for i = 0, nmitx-1 { for j = 0, nmity-1 { for k = 0, mit[i][j].spiketimes.size()-1 { outfile.printf("%d %d %d %10.3f\n",i,j,i*nmity+j,mit[i][j].spiketimes.x[k]) } } } outfile.close() } proc print_gran_raster() { // 1 arg - filename root sprint(filename,"%s.gran.ras",\$s1) outfile.wopen(filename) outfile.printf("# Granule cell raster plot\n") for i = 0, ngranx-1 { for j = 0, ngrany-1 { for k = 0, gran[i][j].spiketimes.size()-1 { outfile.printf("%d %d %d %10.3f\n",i,j,i*ngrany+j,gran[i][j].spiketimes.x[k]) } } } outfile.close() } proc print_isis() { // 3 args - indices of mitral cell plus filename root calc_isis(\$1,\$2,ttrans) sprint(filename,"%s_%d_%d.isi",\$s3,\$1,\$2) outfile.wopen(filename) outfile.printf("# Interspike intervals for mitral cell [%d][%d]",\$1,\$2) isi.printf(outfile,"%10.3f") outfile.close() } proc print_isi_stats() { // 1 arg - filename root sprint(filename,"%s.stats",\$s1) outfile.wopen(filename) outfile.printf("#Interspike interval statistics for mitral cells\n") outfile.printf("# i j n mean median stdev \n") for i = 0, nmitx-1 { for j = 0, nmity-1 { calc_isis(i,j,ttrans) outfile.printf("%3d%3d%4d",i,j,isi.size()) if (isi.size() > 0) { outfile.printf("%8.2f%8.2f",isi.mean(),isi.median()) if (isi.size() > 1) { outfile.printf("%8.2f\n",isi.stdev()) } else { printf("\n") } } else { outfile.printf("\n") } } } outfile.printf("#Interspike interval statistics for granule cells\n") outfile.printf("# i j n mean median stdev \n") for i = 0, ngranx-1 { for j = 0, ngrany-1 { calc_gran_isis(i,j,ttrans) outfile.printf("%3d%3d%4d",i,j,isi.size()) if (isi.size() > 0) { outfile.printf("%8.2f%8.2f",isi.mean(),isi.median()) if (isi.size() > 1) { outfile.printf("%8.2f\n",isi.stdev()) } } else { outfile.printf("\n") } } } outfile.close() } proc print_lags() { local i,j // 3 args - indices of mitral cell + filename root sprint(filename,"%s_%d_%d.lags",\$s3,\$1,\$2) outfile.wopen(filename) outfile.printf("# Lag times for mitral cell [%d][%d]\n",\$1,\$2) for i = 0, nmitx-1 { for j = 0, nmity-1 { calc_lags(\$1,\$2,i,j,ttrans) outfile.printf("[%d,%d]",i,j) lags.printf(outfile,"%10.3f") } } outfile.close() } proc print_si() { // 1 arg - fileroot print "Calculating synchronization indices" sprint(filename,"%s.synch",\$s1) outfile.wopen(filename) outfile.printf("Synchronization index: %10.3f\n",synch_index()) outfile.printf("Phase-locking index: %10.3f\n",phaselock_index()) outfile.close() }