NEURON interfaces to MySQL and the SPUD feature extraction algorithm (Neymotin et al. 2008)

 Download zip file   Auto-launch 
Help downloading and running models
See the readme.txt for information on setting up this interface to a MySQL server from the NEURON simulator. Note the SPUD feature extraction algorithm includes its own readme in the spud directory.
1 . Neymotin S, Uhlrich DJ, Manning KA, Lytton WW (2008) Data mining of time-domain features from neural extracellular field data Applic. of Comput. Intel. in Bioinf. and Biomed.: Current Trends and Open Problems 151:119-140
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s):
Gap Junctions:
Simulation Environment: NEURON;
Model Concept(s): Methods;
Implementer(s): Lytton, William [billl at]; Neymotin, Sam [samn at];
Files displayed below are from the implementation

:SPUD feature extraction algorithm implementation as appearing in:
:Data-mining of time-domain features from neural extracellular field data (2007 - in press)
:S Neymotin, DJ Uhlrich, KA Manning, WW Lytton

  SUFFIX nothing
  : BVBASE is bit vector base number (typically 0 or -1)

  SHM_SPUD=4   : used in spud() for measuring sharpness
  NOV_SPUD=1   : used in spud() to eliminate overlap of spikes
  CREEP_SPUD=0 : used in spud() to allow left/right "creep" to local minima

#include <stdlib.h>
#include <math.h>
#include <values.h> // contains MAXLONG 
#include <sys/time.h> 
extern double* hoc_pgetarg();
extern double hoc_call_func(Symbol*, int narg);
extern FILE* hoc_obj_file_arg(int narg);
extern Object** hoc_objgetarg();
extern void vector_resize();
extern int vector_instance_px();
extern void* vector_arg();
extern double* vector_vec();
extern double hoc_epsilon;
extern double chkarg();
extern void set_seed();
extern int ivoc_list_count(Object*);
extern Object* ivoc_list_item(Object*, int);
extern int hoc_is_double_arg(int narg);
extern char* hoc_object_name(Object*);
char ** hoc_pgargstr();
int list_vector_px();
int list_vector_px2();
int list_vector_px3();
int list_vector_resize();
int ismono1();
static void hxe() { hoc_execerror("",0); }
static void hxf(void *ptr) { free(ptr); hoc_execerror("",0); }

:* src.spud(thresh,dlist,nqslist)
:  dest.spud(src)  -- default thresh=0; returns indices
: look for multiple threshold crossings to define peaks
: creates multiple parallel vectors for an NQS db
: counts peaks pointing upward -- should be all pos
: see eg decnqs.hoc:fudup() for usage

  //** declarations
#define UDSL 50
#define UDNQ 11
#define LOC     nq[0] // loc of peak of spike
#define PEAK  	nq[1] // value at peak (absolute height)
#define WIDTH  	nq[2] // rt flank - lt flanks (? isn't it rt flank - LOC ?)
#define BASE  	nq[3] // height at base
#define HEIGHT  nq[4] // peak - base
#define START  	nq[5] // left flank of spike?
#define SLICES  nq[6] // how many slices found this spike
#define SHARP  	nq[7] // 2nd deriv at peak
#define INDEX  	nq[8] // consecutive numbering of spikes
//        	nq[9] // will use to fill in trace's file name at hoc level
#define NESTED  nq[10] // how many bumps are nested within this one
  //** procedure spud()
static double spud (void* vv) {
  int i, k, m, n, nqsz, nsrc, jj[UDSL], f[UDSL], lc, dsz[UDSL], nqmax, thsz, lc2, done, dbn;
  double *src, *tvec, *th, *dest[UDSL], *nq[UDNQ], *tmp, *dbx, lt, thdist;
  Object *ob, *ob2;
  void *vvd[UDSL], *vvth, *vnq[UDNQ];
  //** read in vectors and verify sizes, etc
  nsrc = vector_instance_px(vv, &src); // trace to analyze
  thsz = vector_arg_px(1, &th);        // vector of thresholds to check
  ob =  *hoc_objgetarg(2);             // storage for values for each threshold
  ob2 = *hoc_objgetarg(3);             // list of NQS vectors for returning values
  tmp = (double *)ecalloc(nsrc, sizeof(double));  // tmp is size of trace
  lc =  ivoc_list_count(ob);
  lc2 = ivoc_list_count(ob2);
  if (lc>UDSL) {printf("spud ERRF mismatch: max slice list:%d %d\n",UDSL,lc); hxf(tmp);}
  if (lc2!=UDNQ){printf("spud ERRB mismatch: NQS sz is %d (%d in list)\n",UDNQ,lc2);hxf(tmp);}
  if (nsrc<lc) {printf("spud ERRC mismatch: %d %d\n",lc,nsrc); hxf(tmp);} // ??
  if (lc!=thsz) {printf("spud ERRA mismatch: %d %d\n",lc,thsz); hxf(tmp);}
  if (!ismono1(th,thsz,-1)) {printf("spud ERRD: not mono dec %g %d\n",th[0],thsz); hxf(tmp);}
  // thdist=(th[thsz-2]-th[thsz-1])/2; // NOT BEING USED: the smallest spike we will accept
  for (k=0;k <lc;k++)  dsz[k] =list_vector_px3(ob , k, &dest[k], &vvd[k]);
  for (k=0;k<lc2;k++) {
    i=list_vector_px3(ob2, k, &nq[k],   &vnq[k]);
    if (k==0) nqmax=i; else if (i!=nqmax) { // all NQ vecs same size
      printf("spud ERRE mismatch: %d %d %d\n",k,i,nqmax); hxf(tmp); }
  //** store crossing points and midpoints in dest[k]
  // dest vectors dest[k] will store crossing points and midpoints at each th[k] slice location
  // as triplets: up/max/down
  for (k=0; k<lc; k++) {   // iterate thru thresholds
    jj[k]=f[k]=0; // jj[k] is ind into dest[k]; f[k] is flag for threshold  crossings
    for (i=0;i<nsrc && src[i]>th[k];i++) {} // start somewhere below this thresh th[k]
    for (; i<nsrc; i++) { // iterate through trace
      if (src[i]>th[k]) { 
        if (f[k]==0) { // ? passing thresh 
          if (jj[k]>=dsz[k]){printf("(%d,%d,%d) :: ",k,jj[k],dsz[k]);
            hoc_execerror("Dest vec too small in spud ", 0); }
          dest[k][jj[k]++] = (i-1) + (th[k]-src[i-1])/(src[i]-src[i-1]); // interpolate
          tmp[k]=-1e9; dest[k][jj[k]]=-1.; // flag in tmp says that a thresh found here
        if (f[k]==1 && src[i]>tmp[k]) { // use tmp[] even more temporarily
          tmp[k]=src[i]; // pick out max
          dest[k][jj[k]] = (double)i; // location of this peak
      } else {          // below thresh 
        if (f[k]==1) {  // just passed going down 
          jj[k]++;      // triplet will be indices of cross-up/peak/cross-down
          dest[k][jj[k]++] = (i-1) + (src[i-1]-th[k])/(src[i-1]-src[i]);
  //** truncate dest vectors to multiples of 3:
  for (k=0;k<lc;k++) vector_resize(vvd[k],(int)(floor((double)jj[k]/3.)*3.));
  for (i=0; i<nsrc; i++) tmp[i]=0.; // clear temp space
  //** go through all the slices to find identical peaks and save widths and locations
  // tmp[] uses triplets centered around a location corresponding to a max loc in the
  // original vector; the widest flanks for each are then on either side of this loc
  for (k=0;k<lc;k++) { // need to go from top to bottom to widen flanks
    for (i=1;i<jj[k];i+=3) { // through centers (peaks)
      m=(int)dest[k][i]; // hash: place center at location
      if (tmp[m-2]<0 || tmp[m-1]<0 || tmp[m+1]<0 || tmp[m+2]<0) continue; // ignore; too crowded
      tmp[m]--;  // count how many slices have found this peak (use negative)
      tmp[m-1]=dest[k][i-1]; tmp[m+1]=dest[k][i+1]; // flanks
  //** 1st (of 2) loops through tmp[] -- pick up flanks
  // step through tmp[] looking for negatives which indicate the slice count and pick up 
  // flanks from these
  for (i=0,k=0; i<nsrc; i++) if (tmp[i]<0.) { // tmp holds neg of count of slices
    if (k>=nqmax) { printf("spud ERRG OOR in NQ db: %d %d\n",k,nqmax); hxf(tmp); }
    LOC[k]=(double)i;  // approx location of the peak of the spike
    WIDTH[k]=tmp[i+1]; // location of right side -- temp storage
    START[k]=tmp[i-1]; // start of spike (left side)
    SLICES[k]=-tmp[i];  // # of slices
  nqsz=k;   // k ends up as size of NQS db
  if (DEBUG_SPUD && ifarg(4)) { dbn=vector_arg_px(4, &dbx); // DEBUG -- save tmp vector
    if (dbn<nsrc) printf("spud ERRH: Insufficient room in debug vec (%d<%d)\n",dbn,nsrc); 
    else for (i=0;i<nsrc;i++) dbx[i]=tmp[i]; 
  //** adjust flanks to handle nested bumps
  // 3 ways to handle spike nested in a spike or elongated base:
  // NB always using same slice for both L and R flanks; NOV_SPUD flag: (no-overlap)
  //   0. nested spike(s) share flanks determined by shared base
  //   1. nested spike(s) have individual bases, 1st and last use flanks from base
  //   2. nested spike(s) have individual bases, base flanks listed separately w/out peak
  // here use 
  // search nq vecs to compare flanks to neighboring centers
  // if flanks overlap the centers on LT or RT side,
  // correct them by going back to original slice loc info (in dest[])
  //*** look at left side -- is this flank to left of center of another bump?
  if (NOV_SPUD) for (i=0;i<nqsz;i++) { // iterate through NQS db
    if ((i-1)>0 && START[i] < LOC[i-1]) { // flank is to left of prior center
      if (DEBUG_SPUD) printf("LT problem %d %g %g<%g\n",i,LOC[i],START[i],LOC[i-1]);
      for (m=lc-1,done=0;m>=0 && !done;m--) { // m:go from bottom (widest) to top
        for (n=1;n<jj[m] && !done;n+=3) {     // n:through centers
          // pick out lowest slice with this peak LOC whose flank is to RT of prior peak
          if (floor(dest[m][n])==LOC[i] && dest[m][n-1]>LOC[i-1]) {
            // ??[i]=START[i]; // temp storage for L end of this overlap
            // replace both left and right flanks at this level -- #1 above
            START[i]=dest[m][n-1]; WIDTH[i]=dest[m][n+1]; done=1; 
    //*** now look at RT side
    if ((i+1)<nqsz && WIDTH[i]>LOC[i+1]) {
      if (DEBUG_SPUD) printf("RT problem %d %g %g>%g\n",i,LOC[i],WIDTH[i],LOC[i+1]);
      for (m=lc-1,done=0;m>=0 && !done;m--) { // m: go from bottom to top
        for (n=1;n<jj[m] && !done;n+=3) {     // n: through centers
          // pick out lowest slice with this peak LOC whose flank is to LT of next peak
          if (floor(dest[m][n])==LOC[i] && dest[m][n+1]<LOC[i+1]) {
            // ??[i]=WIDTH[i]; // end of overlap
            START[i]=dest[m][n-1]; WIDTH[i]=dest[m][n+1]; done=1;

  //make sure left and right sides of bump occur at local minima
  //shouldn't creeping be before NOV_SPUD=1 overlap check???
  //creeping can result only in equal borders btwn two bumps
  //on one side, so it should be ok here...
  if(CREEP_SPUD) for(i=0,k=0;i<nsrc;i++) if(tmp[i]<0.){

    //move left side to local minima
    int idx = (int)START[k];
    while(idx >= 1 && src[idx] >= src[idx-1]) idx--;
    START[k] = idx;

    //move right side to local minima
    idx = (int)WIDTH[k];
    while(idx < nsrc-1 && src[idx] >= src[idx+1]) idx++;
    WIDTH[k] = idx;


  //** 2nd loop through tmp[] used to fill in the rest of NQS
  // needed to split into 2 loops so that could check for overlaps and correct those
  // before filling in the rest of nq
  for (i=0,k=0; i<nsrc; i++) if (tmp[i]<0.) { // tmp holds neg of count of slices
    // calculate a base voltage lt as interpolated value on left side
    BASE[k]=lt;         // base voltage
    PEAK[k]=src[i];     // peak voltage
    WIDTH[k] = WIDTH[k] - START[k]; // width = RT_flank-LT_flank
    HEIGHT[k]=PEAK[k]-BASE[k]; // redund measure -- can eliminate
    // measure of sharpness diff of 1st derivs btwn peak and SHM_SPUD dist from peak
    // to get 2nd deriv would be normalized by 2*SHM_SPUD*tstep
    // ??could take an ave. or max first deriv for certain distance on either side
  int iNumBumps = k;

  //count # of other bumps nested within each bump
    for(i=0; i<iNumBumps; i++){
      NESTED[i] = 0;
      int j = 0;
        if(i!=j && LOC[j] >= START[i] && LOC[j] <= START[i]+WIDTH[i]){
  } else for(i=0;i<iNumBumps;i++) NESTED[i]=0.0;

  //** finish up
  for (i=0;i<lc2;i++) vector_resize(vnq[i], nqsz);
  if (k!=nqsz) { printf("spud ERRI INT ERR: %d %d\n",k,nqsz); hxf(tmp); }
  return jj[0];


:* PROCEDURE install_spud()
PROCEDURE install_spud () {
  if (SPUD_INSTALLED==1) {
    printf("spud.mod installed\n")
  } else {
  install_vector_method("spud", spud);

Neymotin S, Uhlrich DJ, Manning KA, Lytton WW (2008) Data mining of time-domain features from neural extracellular field data Applic. of Comput. Intel. in Bioinf. and Biomed.: Current Trends and Open Problems 151:119-140

References and models cited by this paper

References and models that cite this paper

Adeli H, Zhou Z, Dadmehr N (2003) Analysis of EEG records in an epileptic patient using wavelet transform. J Neurosci Methods 123:69-87 [PubMed]

Ascoli GA, De Schutter E, Kennedy DN (2003) An information science infrastructure for neuroscience. Neuroinformatics 1:1-2 [Journal] [PubMed]

Carnevale NT, Hines ML (2006) The NEURON Book

Chamberlin DD,Boyce RF (1974) SEQUEL: A structured English query language, International Conference on Management of Data Proceedings of the 1974 ACM SIGFIDET (now SIGMOD) workshop on Data description, access and control, Ann Arbor, Michigan :249-264

Cooley JW, Connor CE (1965) An algorithm for the machine calculation of complex Fourier series Math Comput 19:297-301

Dunn JC (1973) A fuzzy relative of the ISODATA process and its use in detecting compact Well-Separated clusters Journal of Cybernetics 3:32-57 [Journal]

Golarai G, Cavazos JE, Sutula TP (1992) Activation of the dentate gyrus by pentylenetetrazol evoked seizures induces mossy fiber synaptic reorganization. Brain Res 593:257-64 [PubMed]

Hartigan JA, Wong MA (1979) Algorithm as 136: a k-means clustering algorithm J R Stat Soc Ser C 28:100-108

Hines ML, Carnevale NT (1997) The NEURON simulation environment. Neural Comput 9:1179-209 [PubMed]

Hines ML, Carnevale NT (2000) Expanding NEURON's repertoire of mechanisms with NMODL. Neural Comput 12:995-1007 [PubMed]

Johnson SG,Frigo M (2007) A modified split-radix FFT with fewer arithmetic operations IEEE Transactions on Signal Processing 55:111-119

Lytton WW (2006) Neural Query System: Data-mining from within the NEURON simulator. Neuroinformatics 4:163-76 [Journal] [PubMed]

   Neural Query System NQS Data-Mining From Within the NEURON Simulator (Lytton 2006) [Model]

Lytton WW, Omurtag A (2007) Tonic-clonic transitions in computer simulation. J Clin Neurophysiol 24:175-81 [PubMed]

   Tonic-clonic transitions in a seizure simulation (Lytton and Omurtag 2007) [Model]

Macqueen JB (1967) Some methods for classification and analysis of multivariate observations Proceedings Of 5th Berkeley Symposium On Mathematical Statistics And Probability 1:281-297

Sweldens W (1997) The lifting scheme: A construction of second generation wavelets SIAM Journal on Mathematical Analysis 29:511-546

Uhlrich DJ, Manning KA, O'Laughlin ML, Lytton WW (2005) Photic-induced sensitization: acquisition of an augmenting spike-wave response in the adult rat through repeated strobe exposure. J Neurophysiol 94:3925-37 [PubMed]

Neymotin SA, Hilscher MM, Moulin TC, Skolnick Y, Lazarewicz MT, Lytton WW (2013) Ih Tunes Theta/Gamma Oscillations and Cross-Frequency Coupling In an In Silico CA3 Model PLoS ONE 8(10):e76285 [Journal] [PubMed]

   Ih tunes oscillations in an In Silico CA3 model (Neymotin et al. 2013) [Model]

Neymotin SA, Lee H, Park E, Fenton AA, Lytton WW (2011) Emergence of physiological oscillation frequencies in a computer model of neocortex. Front Comput Neurosci 5:19-75 [Journal] [PubMed]

   Emergence of physiological oscillation frequencies in neocortex simulations (Neymotin et al. 2011) [Model]

(18 refs)