N-Neuron Network Coherence Code

Posted by Will on June 26th 2006 to Neuroscience

Here's some code I wrote to take spiketimes and find the average coherence between each pair of neurons. Files must be named xvoltage_raster.out, and each line should start with the time of spike for the xth neuron. The directory containing the raster files should be in this format: <neurons>_<gsyn>_<iapp>_<tausyn>_<heterogeneity>. An example would be: /50_0.35_3.0_1.0_0.038/. The method is from White et al., 1998, and the epsilon/heterogeneity function comes from Skinner et al., 2005. This is a rough draft, but if I don't post it now I'll probably forget later. Click through for the code.

Warning: This code determines the window size by the fastest neuron in the network, not by the fastest of each respective pair. The literature isn't too clear, and I don't have experience with cross-correlation outside of this application. If you want the 'correct' code, contact me. Otherwise the concept is similar except for the computing of each pair separately.

PLAIN TEXT
C:
  1. #include <stdio.h>
  2. #include <pthread.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include <unistd.h>
  6.  
  7. double eps2het(double iapp, double epsilon)
  8. {
  9.     double het, iapp1, iapp2, intfreq1, intfreq2;
  10.  
  11.     iapp1 = iapp - epsilon;
  12.     iapp2 = iapp + epsilon;
  13.  
  14.     intfreq1 = (-0.112*pow(iapp1,4)) + (1.88*pow(iapp1,3)) - (13.1*pow(iapp1,2)) + (70.5*iapp1) + 0.258;
  15.     intfreq2 = (-0.112*pow(iapp2,4)) + (1.88*pow(iapp2,3)) - (13.1*pow(iapp2,2)) + (70.5*iapp2) + 0.258;
  16.  
  17.     het = ((intfreq2 - intfreq1)/intfreq2);
  18.  
  19.     return het;
  20. }
  21.  
  22. double het2eps(double iapp, double hetero)
  23. {
  24.     double testeps = 0.1;
  25.     double step = 0.000001;
  26.     while ( 1 )
  27.     {
  28.         double newhet = eps2het(iapp, testeps);
  29.         double hetdiff = newhet - hetero;
  30.         double hetdiffabs = fabs(hetdiff);
  31.         //printf("%lf\n",hetdiffabs);
  32.  
  33.         if (hetdiffabs <step)
  34.         {
  35.             return testeps;
  36.         }
  37.         else if (hetdiff> 0)
  38.             testeps = testeps - step;
  39.         else if (hetdiff <0)
  40.             testeps = testeps + step;
  41.     }
  42. }
  43.  
  44. double unitsquarewavearea(short *waveform, int length, double binsize)
  45. {
  46.     double area = 0;
  47.  
  48.     int i;
  49.     for (i = 0; i <length; i++)
  50.     {
  51.         if (waveform[i] == 1)
  52.         {
  53.             area = area + binsize;
  54.         }
  55.     }
  56.  
  57.     return area;
  58. }
  59.  
  60. void setrelativedir(char *reldir)
  61. {
  62.     int i = 0;
  63.     int j = 0;
  64.  
  65.     char cdir[200];
  66.     getcwd(cdir,sizeof(char) * 200);
  67.  
  68.     for (i = 0; cdir[i]> 0; i++)
  69.     {
  70.         if (cdir[i] == '/')
  71.         {
  72.             reldir[j] = 0;
  73.             j = 0;
  74.         }
  75.         else
  76.         {
  77.             reldir[j] = cdir[i];
  78.             j++;
  79.         }
  80.  
  81.         reldir[j] = 0;
  82.     }
  83. }
  84.  
  85. void setvaluefromdir(char *dir, int *neurons, double *gsyn, double *iapp, double *tausyn, double *het)
  86. {
  87.     int i = 0;
  88.     int j = 0;
  89.     int par = 0;
  90.     int done = 0;
  91.  
  92.     double parval;
  93.     char *strparam = malloc(200);;
  94.  
  95.     for (i = 0; done != 1 ; i++)
  96.     {
  97.         if (dir[i] == '\0' && i> 0 && par> 5)
  98.         {
  99.             done = 1;
  100.         }
  101.  
  102.         if (dir[i] == '_' || (dir[i] == '\0' && i> 0))
  103.         {
  104.             char *endptr;
  105.             strparam[j] = 0;
  106.  
  107.             parval = strtod(strparam, &endptr);
  108.  
  109.             j = 0;
  110.  
  111.             switch (par)
  112.             {
  113.                 case 0:
  114.                     *neurons = (int)parval;
  115.                     break;
  116.                 case 1:
  117.                     *gsyn = parval;
  118.                     break;
  119.                 case 2:
  120.                     *iapp = parval;
  121.                     break;
  122.                 case 3:
  123.                     *tausyn = parval;
  124.                     break;
  125.                 case 4:
  126.                     *het = parval;
  127.                     break;
  128.                 default:
  129.                     break;
  130.             }
  131.  
  132.             par++;
  133.         }
  134.         else
  135.         {
  136.             strparam[j] = dir[i];
  137.             j++;
  138.         }
  139.     }
  140.  
  141.     free(strparam);
  142. }
  143.  
  144. int main(int argc, char *argv[])
  145. {
  146.     int neurons;
  147.     int *pneurons;
  148.     pneurons = &neurons;
  149.  
  150.     int maxspikes, totalbins, leftshift;
  151.     int i, j, k;
  152.  
  153.     double epsilon, stepsize, windowsize, starttime, endtime, timelength;
  154.     double itotal, fastfreq, periodms, window, halfwindow;
  155.  
  156.     double gsyn, iapp, tausyn, het;
  157.     double *pgsyn = &gsyn, *piapp = &iapp, *ptausyn = &tausyn, *phet = &het;
  158.  
  159.     char reldir[200];
  160.     setrelativedir(reldir);
  161.  
  162.     //printf("reldir: %s\n", reldir);
  163.  
  164.     setvaluefromdir(reldir, pneurons, pgsyn, piapp, ptausyn, phet);
  165.  
  166.     epsilon = het2eps(iapp, het);
  167.  
  168.     printf("neurons: %d\niapp: %lf\nhet: %lf\nepsilon: %lf\n\n",neurons,iapp,het,epsilon);
  169.  
  170.     stepsize = 0.01;
  171.     windowsize = 0.2;
  172.  
  173.     starttime = 7000.0;
  174.     endtime = 8000.0;
  175.  
  176.     timelength = endtime - starttime;
  177.  
  178.     totalbins = (int)ceil(timelength / stepsize);
  179.     //printf("bins: %d\n", totalbins);
  180.  
  181.     itotal = iapp + epsilon;
  182.  
  183.     fastfreq = -0.112*pow(itotal,4) + 1.88*pow(itotal,3) - 13.1*pow(itotal,2) + 70.5*itotal + 0.258;
  184.  
  185.     periodms = 1000 / fastfreq;
  186.     //printf("periodms: %lf\n", periodms);
  187.  
  188.     window = windowsize * periodms;
  189.     halfwindow = 0.5 * window;
  190.     leftshift = (int)ceil(halfwindow/stepsize);
  191.     //printf("leftshift: %d\n", leftshift);
  192.  
  193.     FILE *vraster[neurons];
  194.  
  195.     maxspikes = (int)ceil((timelength/periodms) * 2);
  196.     //printf("%d\n", maxspikes);
  197.     double spiketrain[neurons][maxspikes];
  198.     double area[neurons];
  199.  
  200.     short *waveform = malloc(sizeof(short) * neurons * totalbins);
  201.  
  202.     for (i = 0; i <neurons; i++)
  203.     {
  204.         char filename[42];
  205.         sprintf(filename, "%dvoltage_raster.out", i+1);
  206.         //printf("%d %s\n", i+1, filename);
  207.         vraster[i] = fopen(filename, "r");
  208.  
  209.         j = 0;
  210.         double tempdouble;
  211.         while (fscanf(vraster[i], "%lf   %d", &tempdouble, &k) == 2)
  212.         {
  213.             if (tempdouble> starttime && tempdouble <= endtime)
  214.             {
  215.                 spiketrain[i][j] = tempdouble;
  216.                 j++;
  217.             }
  218.         }
  219.  
  220.         fclose(vraster[i]);
  221.  
  222.         for (j = 0; j <totalbins; j++)
  223.         {
  224.             double spiketest = starttime + j*stepsize;
  225.             waveform[(i*totalbins)+j] = 0;
  226.  
  227.             int m;
  228.             for (m = 0; m <maxspikes; m++)
  229.             {
  230.                 if (spiketrain[i][m]> 0)
  231.                 {
  232.                     double timediff = fabs(spiketest - spiketrain[i][m]);
  233.  
  234.                     if (timediff <0.001)
  235.                     {
  236.                         int n;
  237.                         for (n = 0; n <(leftshift*2); n++)
  238.                         {
  239.                             int postshift = (j - leftshift) + n;
  240.                             if (postshift> 0 && postshift <totalbins)
  241.                             {
  242.                                 waveform[(i*totalbins)+postshift] = 1;
  243.                             }
  244.                         }
  245.                         j = j + leftshift - 1;
  246.                     }
  247.                 }
  248.             }
  249.         }
  250.  
  251.         area[i] = unitsquarewavearea(waveform + i*totalbins, totalbins, stepsize);
  252.     }
  253.  
  254.     short *product = malloc(sizeof(short) * totalbins);
  255.     double rhosum = 0;
  256.     int perms = 0;
  257.     int x, y, z;
  258.     for (x = (neurons - 2); x>= 0; x--)
  259.     {
  260.         for (y = (neurons - 1); y> x; y--)
  261.         {
  262.             for (z = 0; z <totalbins; z++)
  263.             {
  264.                 product[z] = waveform[(x*totalbins)+z] * waveform[(y*totalbins)+z];
  265.             }
  266.  
  267.             double productarea = unitsquarewavearea(product, totalbins, stepsize);
  268.             double rho = productarea / (sqrt(area[x]*area[y]));
  269.  
  270.             rhosum = rhosum + rho;
  271.  
  272.             perms++;
  273.             //printf("area_%d: %lf area_%d: %lf rho: %lf\n", x, area[x], y, area[y], rho);
  274.         }
  275.     }
  276.  
  277.     double avgrho = rhosum / perms;
  278.  
  279.     printf("average rho: %lf\n", avgrho);
  280.  
  281.     free(waveform);
  282.     free(product);
  283.  
  284.     return 0;
  285. }

4 Responses to “N-Neuron Network Coherence Code”

  1. ryan
    July 4th, 2006 | 10:56 am

    http://www.brainbasedbusiness.com/2006/07/can_science_really_blueprint_t.html

  2. Todd
    July 27th, 2006 | 7:54 pm

    I hate average rhos. I tend to go for the exceptional rhos, myself.

    I JUST POSTED GARBAGE IN YOUR DISCUSSION FORUM!!

  3. Dillon
    September 5th, 2006 | 6:32 am

    thats newbie code.

  4. Kevs
    August 5th, 2007 | 1:38 pm

    Yowtch! (O_o)~

Leave a reply