Thalamocortical augmenting response (Bazhenov et al 1998)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:37819
In the cortical model, augmenting responses were more powerful in the "input" layer compared with those in the "output" layer. Cortical stimulation of the network model produced augmenting responses in cortical neurons in distant cortical areas through corticothalamocortical loops and low-threshold intrathalamic augmentation. ... The predictions of the model were compared with in vivo recordings from neurons in cortical area 4 and thalamic ventrolateral nucleus of anesthetized cats. The known intrinsic properties of thalamic cells and thalamocortical interconnections can account for the basic properties of cortical augmenting responses. See reference for details. NEURON implementation note: cortical SU cells are getting slightly too little stimulation - reason unknown.
Reference:
1 . Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ (1998) Computational models of thalamocortical augmenting responses. J Neurosci 18:6444-65 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Thalamus;
Cell Type(s): Thalamus geniculate nucleus (lateral) principal neuron; Thalamus reticular nucleus cell; Neocortex V1 pyramidal corticothalamic L6 cell;
Channel(s): I Na,t; I T low threshold; I A; I K,Ca;
Gap Junctions:
Receptor(s): GabaA; GabaB; AMPA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Synchronization; Synaptic Integration;
Implementer(s): Lytton, William [billl at neurosim.downstate.edu];
Search NeuronDB for information about:  Thalamus geniculate nucleus (lateral) principal neuron; Thalamus reticular nucleus cell; Neocortex V1 pyramidal corticothalamic L6 cell; GabaA; GabaB; AMPA; I Na,t; I T low threshold; I A; I K,Ca; Gaba; Glutamate;
: $Id: vecst.mod,v 1.144 2004/07/10 19:07:01 billl Exp $
 
:* COMMENT
COMMENT
randwd   randomly chooses n bits to set to 1
hamming  v.hamming(v1) is hamming distance between 2 vecs
flipbits v.flipbits(scratch,num) flips num rand chosen bits
flipbalbits v.flipbalbits(scratch,num) balanced flipping
vpr      v.vpr prints out vector as 1 (x[i]>0) or 0 (x[i]<=0)
thresh   turns analog vec to BVBASE,1 vec separating at thresh (scalar or vec)
triplet  return location of a triplet in a vector
onoff    turns state vec on or off depending on values in other vecs
bpeval   used by backprop: vo=outp*(1-outp)*del
w        like where but sets chosen nums // modified from .wh in 1.23
xing     determines where vector crosses in a positive direction (destructive)
xzero    count 0 crossings by putting in a 1 or -1
indset(ind,val) sets spots indexed by ind to val
ismono([arg]) checks if a vector is monotonically increaseing (-1 decreasing)
count(num) count the number of nums in vec
scr.fewind(ind,vec1,vec2[,vec3,vec4]) // uses ind as index into other vecs
ind.findx(vecA1,vecA2,...,vecB1,vecB2,...) // uses ind as index into vecA's to load in vecB's
ind.sindx(vecA1,vecA2,...,vecB1,vecB2,...) // replace ind elements in vecA's with vecB's values
ind.sindv(vecA1,vecA2,...,valvec) // replace ind elements in vecA's with vals
scr.nind(ind,vec1,vec2[,vec3,vec4]) // uses ind to elim elems in other vecs
ind.keyind(key,vec1,vec2[,vec3,vec4]) // pick out bzw. values from vectors
ind.slct(key,args,vec1,vec2[,vec3,vec4]) // pick out bzw. values from vectors
vdest.intrp(flag) // interpolate numbers replacing numbers given as flag
v.insct(v1,v2)   // return v1 intersect v2
vdest.cull(vsrc,key)   // remove values found in key
vdest.redundout(vsrc[,INDFLAG])  // remove repeat values, with flag return indices
v.cvlv(v1,v2)   // convolve v1 with v2
fac      not vec related - returns factorial
logfac   not vec related - returns log factorial
v2d      copy into a double block
d2v      copy out of a double block NB: same as vec.from_double()
smgs     rewrite of sumgauss form ivovect.cpp
iwr      write integers
ird      read integers
ENDCOMMENT

NEURON {
  SUFFIX nothing
  GLOBAL BVBASE, RES, VECST_INSTALLED  : bit vector base number (typically 0 or -1)
}

PARAMETER {
  BVBASE = -1.
  VECST_INSTALLED=0
  : misc
  ERR=-1.3479e121

  : 0 args
  ALL=-1.3479e120
  NEG=-1.3478e120
  POS=-1.3477e120
  CHK=-1.3476e120
  NOZ=-1.3475e120
  : 1 arg
  GTH=-1.3474e120
  GTE=-1.3473e120
  LTH=-1.3472e120
  LTE=-1.3471e120
  EQU=-1.3470e120
  EQV=-1.3469e120
  NEQ=-1.3468e120
  SEQ=-1.3467e120
  RXP=-1.3466e120
  : 2 args
  IBE=-1.3465e120
  EBI=-1.3464e120
  IBI=-1.3463e120
  EBE=-1.3462e120
}

ASSIGNED { RES }

VERBATIM
#include <stdlib.h>
#include <math.h>
/* #include <values.h> /* contains MAXLONG */
#include <limits.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 void vector_resize();
extern int vector_instance_px();
extern void* vector_arg();
extern double hoc_epsilon;
ENDVERBATIM
 
:* v1.slope(num) does a linear regression to find the slope, assuming num=timestep of vector

VERBATIM
static double slope(void* vv) {
	int i, n;
	double *x, *y;
        double timestep, sigxy, sigx, sigy, sigx2;
	/* how to get the instance data */
	n = vector_instance_px(vv, &y);

        if(ifarg(1)) { 
          timestep = *getarg(1); 
        } else { printf("You must supply a timestep\n"); return 0; }

        sigxy= sigx= sigy= sigx2=0; // initialize these

        x = (double *) malloc(sizeof(double)*n);
        for(i=0; i<n; i++) {
          x[i] = timestep*i;
          sigxy += x[i] * y[i];
          sigx  += x[i];
          sigy  += y[i];
          sigx2 += x[i]*x[i];
        }
        return (n*sigxy - sigx*sigy)/(n*sigx2 - sigx*sigx);
}
ENDVERBATIM
 
:* v1.vslope(v2) does a linear regression, using v2 as the x-coords

VERBATIM
static double vslope(void* vv) {
	int i, n;
	double *x, *y;
        double timestep, sigxy, sigx, sigy, sigx2;
	/* how to get the instance data */
	n = vector_instance_px(vv, &y);

        if(ifarg(1)) {
          if(vector_arg_px(1, &x) != n ) {
            hoc_execerror("Vector size doesn't match.", 0); 
          }
          sigxy= sigx= sigy= sigx2=0; // initialize these

          for(i=0; i<n; i++) {
            sigxy += x[i] * y[i];
            sigx  += x[i];
            sigy  += y[i];
            sigx2 += x[i]*x[i];
          }
        }         
        return (n*sigxy - sigx*sigy)/(n*sigx2 - sigx*sigx);
}
ENDVERBATIM
 
:* v1.stats(num) does a linear regression, assuming num=timestep of vector

VERBATIM
static double stats(void* vv) {
	int i, n;
	double *x, *y;
        double timestep, sigxy, sigx, sigy, sigx2, sigy2;
        double r, m, b;
	/* how to get the instance data */
	n = vector_instance_px(vv, &y);

        if(ifarg(1)) { 
          timestep = *getarg(1); 
        } else { printf("You must supply a timestep\n"); return 0; }

        sigxy= sigx= sigy= sigx2=sigy2= 0; // initialize these

        x = (double *) malloc(sizeof(double)*n);
        for(i=0; i<n; i++) {
          x[i] = timestep*i;
          sigxy += x[i] * y[i];
          sigx  += x[i];
          sigy  += y[i];
          sigx2 += x[i]*x[i];
          sigy2 += y[i]*y[i];
        }
        m = (n*sigxy - sigx*sigy)/(n*sigx2 - sigx*sigx);
        b = (sigy*sigx2 - sigx*sigxy)/(n*sigx2 - sigx*sigx);
        r = (n*sigxy - sigx*sigy)/(sqrt(n*sigx2-sigx*sigx) * sqrt(n*sigy2-sigy*sigy));

        printf("Examined %d data points\n", n);
        printf("slope     = %f\n", m);
        printf("intercept = %f\n", b);
        printf("R         = %f\n", r);
        printf("R-squared = %f\n", r*r);
        return 1;
}
ENDVERBATIM
 
:* v1.vstats(v2) does a linear regression, using v2 as the x-coords

VERBATIM
static double vstats(void* vv) {
	int i, n;
	double *x, *y;
        double timestep, sigxy, sigx, sigy, sigx2, sigy2;
        double r, m, b;
	/* how to get the instance data */
	n = vector_instance_px(vv, &y);

        if(ifarg(1)) {
          if(vector_arg_px(1, &x) != n ) {
            hoc_execerror("Vector size doesn't match.", 0); 
          }
          sigxy= sigx= sigy= sigx2=sigy2=0; // initialize these

          for(i=0; i<n; i++) {
            sigxy += x[i] * y[i];
            sigx  += x[i];
            sigy  += y[i];
            sigx2 += x[i]*x[i];
            sigy2 += y[i]*y[i];
          }
          m = (n*sigxy - sigx*sigy)/(n*sigx2 - sigx*sigx);
          b = (sigy*sigx2 - sigx*sigxy)/(n*sigx2 - sigx*sigx);
          r = (n*sigxy - sigx*sigy)/(sqrt(n*sigx2-sigx*sigx) * sqrt(n*sigy2-sigy*sigy));

          printf("Examined %d data points\n", n);
          printf("slope     = %f\n", m);
          printf("intercept = %f\n", b);
          printf("R         = %f\n", r);
          printf("R-squared = %f\n", r*r);
          return 1;
        } else {
          printf("You must supply an x vector!\n");
          return 0;
        }
}
ENDVERBATIM
 
:* v1.randwd(num[,v2]) will randomly flip num bits from BVBASE to 1
: does v1.fill(BVBASE); optionally fill v2 with the indices
VERBATIM
static double randwd(void* vv) {
	int i, ii, jj, nx, ny, flip, flag;
	double* x, *y;
	/* how to get the instance data */
	nx = vector_instance_px(vv, &x);
        flip = (int) *getarg(1);
        if (ifarg(2)) { /* write a diff vector to z */
          flag = 1; ny = vector_arg_px(2, &y);
          if (ny!=flip) { hoc_execerror("Opt vector must be size for # of flips", 0); }
        } else { flag = 0; }
        if (flip>=nx) { hoc_execerror("# of flips exceeds (or ==) vector size", 0); }
	for (i=0; i < nx; i++) { x[i] = BVBASE; }
	for (i=0,jj=0; i < flip; i++) { /* flip these bits */
	  ii = (int) ((nx+1)*drand48());
	  if (x[ii]==BVBASE) {
	    x[ii] = 1.; 
            if (flag) { y[jj] = ii; jj++; }
	  } else {
	    i--;
	  }
	}
	return flip;
}
ENDVERBATIM
 
:* v1.hamming(v2[,v3]) compares v1 and v2 for matches, v3 gives diff vector
VERBATIM
static double hamming(void* vv) {
	int i, nx, ny, nz, prflag;
	double* x, *y, *z,sum;
	sum = 0.;
	nx = vector_instance_px(vv, &x);
	ny = vector_arg_px(1, &y);
        if (ifarg(2)) { /* write a diff vector to z */
          prflag = 1; nz = vector_arg_px(2, &z);
        } else { prflag = 0; }
	if (nx!=ny || (prflag && nx!=nz)) {
	  hoc_execerror("Vectors must be same size", 0);
	}
	for (i=0; i < nx; ++i) {
	  if (x[i] != y[i]) { sum++; 
            if (prflag) { z[i] = 1.; }
          } else if (prflag) { z[i] = 0.; }
        }
	return sum;
}
ENDVERBATIM
 
:* v1.indset(ind,x[,y]) sets indexed values to x and other values to optional y
VERBATIM
static double indset(void* vv) {
	int i, nx, ny, nz, flag;
	double* x, *y, *z, val, val2 ;
	nx = vector_instance_px(vv, &x);
	ny = vector_arg_px(1, &y);
        if (hoc_is_object_arg(2)) { 
          flag=1;
          nz = vector_arg_px(2, &z); 
          if (ny!=nz) { hoc_execerror("v.indset: Vector sizes don't match.", 0); }
        } else { flag=0; val=*getarg(2); }
        if (ifarg(3)) { 
          val2 = *getarg(3); 
          for (i=0; i<nx; i++) { x[i]=val2; }
        }
	for (i=0; i<ny; i++) {
	  if (y[i] > nx) { hoc_execerror("v.indset: Index exceeds vector size", 0); }
          if (flag) x[(int)y[i]]=z[i]; else x[(int)y[i]]=val;
        }
	return i;
}
ENDVERBATIM
 
VERBATIM
/* Maintain parallel int vector to avoid slowness of repeated casts */
static int scrsz=0;
static int *scr;
static double dcr[100]; // scratch area for doubles
ENDVERBATIM

:* tmp.fewind(ind,vec1,vec2[,vec3,vec4,...])
: picks out numbers from multiple vectors using index ind
VERBATIM
static double fewind(void* vv) {
	int i, j, nx, ni, nv[10], num;
	double *x, *ind, *vvo[10];
	nx = vector_instance_px(vv, &x);
        for (i=1;ifarg(i);i++);
        i--; /* counted 1 arg too many */
	if (i>11) hoc_execerror("ERR: fewind can only handle 10 vectors", 0);
	num = i-1; /* number of vectors to be picked apart */
        for (i=0;i<num;i++) { 
          nv[i] = vector_arg_px(i+2, &vvo[i]);
          if (nx!=nv[i]) { printf("fewind ERR %d %d %d\n",i,nx,nv[i]);
            hoc_execerror("Vectors must all be same size: ", 0); }
        }
        ni = vector_arg_px(1, &ind);
        if (ni>scrsz) { 
          if (scrsz>0) { free(scr); scr=(int *)NULL; }
          scrsz=ni+10000;
          scr=(int *)ecalloc(scrsz, sizeof(int));
        }
        for (i=0;i<ni;i++) { scr[i]=(int)ind[i]; /* copy into integer array */
          if (scr[i]>nx || scr[i]<0) { printf("fewind ERR1 %d %d\n",scr[i],nx);
            hoc_execerror("Index vector out-of-bounds", 0); }
        }
        for (j=0;j<num;j++) {
          for (i=0;i<ni;i++) x[i]=vvo[j][scr[i]];    
          for (i=0;i<ni;i++) vvo[j][i]=x[i];   
          vv=vector_arg(j+2); vector_resize(vv, ni);
        }
	return ni;
}
ENDVERBATIM


:* ind.findx(vecA1,vecA2,...,vecB1,vecB2,...)) 
: uses ind as index into vecA's to load in vecB's (for select); a nondestructive fewind()
VERBATIM
static double findx(void* vv) {
	int i, j, ni, nx, av[11], bv[11], num;
	double *ind, *avo[11], *bvo[11];
	ni = vector_instance_px(vv, &ind);
        for (i=1;ifarg(i);i++);
        i--; /* counted 1 arg too many */
        if (i>22) hoc_execerror("findx ****ERRA****: can only handle 11 vectors", 0);
        if ((i+1)/2 != (i/2)) { printf("\t** %d vecs\n",i);
          hoc_execerror("findx ****ERRB**** requires even # of vecs", 0); }
        num = i/2; /* number of vectors to be picked apart */
        for (i=0;i<num;i++) { 
          av[i]=vector_arg_px(i+1, &avo[i]); // source vectors
          if (av[0]!=av[i]) { printf("findx ****ERRC**** %d %d %d\n",i,av[0],av[i]);
            hoc_execerror("Src vectors must all be same size: ", 0); }
        }
        nx=av[0]; // size of source vecs
        for (i=0;i<num;i++) { 
          bv[i]=vector_arg_px(num+i+1, &bvo[i]); /* dest vectors */
          vv=vector_arg(num+i+1);
          if (vector_buffer_size(vv)<ni) { 
            printf("findx ****ERRD**** arg#%d need:%d sz:%d\n",num+i+1,ni,vector_buffer_size(vv));
            hoc_execerror("Destination vector with insufficient size: ", 0); }
        }
        if (ni>scrsz) { 
          if (scrsz>0) { free(scr); scr=(int *)NULL; }
          scrsz=ni+10000;
          scr=(int *)ecalloc(scrsz, sizeof(int));
        }
        for (i=0;i<ni;i++) { scr[i]=(int)ind[i]; // copy into integer array
          if (scr[i]>nx || scr[i]<0) { printf("findx ****ERRE**** **** IND:%d SZ:%d\n",scr[i],nx);
            hoc_execerror("Index vector out-of-bounds", 0); }
        }
        for (j=0,i=0;j<num;j++) { // j goes through vectors
          vv=vector_arg(j+num+1); 
          vector_resize(vv, ni);
          for (i=0;i<ni;i++) bvo[j][i]=avo[j][scr[i]];    
        }
	return ni;
}
ENDVERBATIM

:* ind.sindx(vecA1,vecA2,...,vecB1,vecB2,...)) 
: uses ind as index into vecA's to replace with elements from vecB's
VERBATIM
static double sindx(void* vv) {
	int i, j, ni, nx, av[11], bv[11], num;
	double *ind, *avo[11], *bvo[11];
	ni = vector_instance_px(vv, &ind);
        for (i=1;ifarg(i);i++);
        i--; /* counted 1 arg too many */
        if (i>22) hoc_execerror("sindx ****ERRA****: can only handle 11 vectors", 0);
        if ((i+1)/2 != (i/2)) { printf("\t** %d vecs\n",i);
          hoc_execerror("sindx ****ERRB**** requires even # of vecs", 0); }
        num = i/2; /* number of vectors to be picked apart */
        for (i=0;i<num;i++) { 
          av[i]=vector_arg_px(i+1, &avo[i]); // dest vectors
          if (av[0]!=av[i]) { printf("sindx ****ERRC**** %d %d %d\n",i,av[0],av[i]);
            hoc_execerror("Dest. vectors must all be same size: ", 0); }
        }
        nx=av[0]; // size of dest vecs
        for (i=0;i<num;i++) { 
          bv[i]=vector_arg_px(num+i+1, &bvo[i]); /* source vectors */
          if (bv[i]!=ni) { 
            printf("sindx ****ERRD**** arg#%d does note match ind length %d vs %d\n",num+i+1,ni,bv[i]);
            hoc_execerror("Source vector with insufficient size: ", 0); }
        }
        if (ni>scrsz) { 
          if (scrsz>0) { free(scr); scr=(int *)NULL; }
          scrsz=ni+10000;
          scr=(int *)ecalloc(scrsz, sizeof(int));
        }
        for (i=0;i<ni;i++) { scr[i]=(int)ind[i]; // copy into integer array
          if (scr[i]>nx || scr[i]<0) { 
            printf("sindx ****ERRE**** IND:%d SZ:%d\n",scr[i],nx);
            hoc_execerror("Index vector out-of-bounds", 0); }
        }
        for (j=0,i=0;j<num;j++) for (i=0;i<ni;i++) avo[j][scr[i]]=bvo[j][i];
	return ni;
}

ENDVERBATIM
:* ind.sindv(vecA1,vecA2,...,valvec)
: uses ind as index into vecA's to replace with values in valvec
VERBATIM
static double sindv(void* vv) {
	int i, j, ni, nx, av[11], bv, num;
	double *ind, *avo[11], *bvo;
	ni = vector_instance_px(vv, &ind);
        for (i=1;ifarg(i);i++);
        num = i-1; /* total number of arguments */
        if (num>12) hoc_execerror("sindv ****ERRA****: can only handle 11 vectors", 0);
        for (i=0;i<num-1;i++) { 
          av[i]=vector_arg_px(i+1, &avo[i]); // dest vectors
          if (av[0]!=av[i]) { printf("sindv ****ERRC**** %d %d %d\n",i,av[0],av[i]);
            hoc_execerror("Dest. vectors must all be same size: ", 0); }
        }
        nx=av[0]; // size of source vecs
        bv=vector_arg_px(num, &bvo); // source vector
        if (bv!=num-1) { 
          printf("sindv ****ERRD**** arg#%d does note match length %d vs %d\n",num+1,num,bv);
          hoc_execerror("Source vector with insufficient size: ", 0); }
        if (ni>scrsz) { 
          if (scrsz>0) { free(scr); scr=(int *)NULL; }
          scrsz=ni+10000;
          scr=(int *)ecalloc(scrsz, sizeof(int));
        }
        for (i=0;i<ni;i++) { scr[i]=(int)ind[i]; // copy into integer array
          if (scr[i]>nx || scr[i]<0) { 
            printf("sindv ****ERRE**** IND:%d SZ:%d\n",scr[i],nx);
            hoc_execerror("Index vector out-of-bounds", 0); }
        }
        for (j=0,i=0;j<num-1;j++) for (i=0;i<ni;i++) avo[j][scr[i]]=bvo[j];
	return ni;
}
ENDVERBATIM

:* ind.slct(key,args,vec1,vec2[,vec3,vec4,...])
: picks out indices of numbers in key from multiple vectors
VERBATIM
static double slct(void* vv) {
  int i, j, k, m, n, ni, nk, na, nv[10], num, fl;
  double *ind, *key, *arg, *vvo[10];
  ni = vector_instance_px(vv, &ind); // vv is ind
  for (i=0;ifarg(i);i++); i--; // drop back by one to get numarg()
  if (i>12) hoc_execerror("ERR: vecst::slct can only handle 10 vectors", 0);
  num = i-2; /* number of vectors to be picked apart */
  for (i=0;i<num;i++) { 
    nv[i] = vector_arg_px(i+3, &vvo[i]);
    if (ni!=nv[i]) { printf("vecst::slct ERR %d %d %d\n",i,ni,nv[i]);
    hoc_execerror("index and searched vectors must all be same size: ", 0); }
  }
  nk = vector_arg_px(1, &key);
  na = vector_arg_px(2, &arg);
  if (2*nk!=na) { printf("vecst::slct ERR3 %d %d\n",nk,na); 
    hoc_execerror("Arg vector must be double key length",0); }
  for (i=0,n=0;i<nk;i++) if (key[i]==EQV) n++; // special case takes 2 vec args
  if (nk+n!=num) { 
    printf("vecst::slct ERR2 %d(keys)+%d(EQV)!=%d(vecs)\n",nk,n,num); 
    hoc_execerror("Key length must be number of vecs + num of EQV",0); }
  for (j=0,k=0,m=0;j<ni;j++) { // j steps through elements of vectors
    for (i=0,m=0,n=0,fl=1;i<num;i++,n++,m+=2) { // i steps thru key, m thru args
      if (key[n]==ALL) continue; // OK - do nothing
      if (key[n]==NOZ)        { if (vvo[i][j]==0.) {fl=0; break;} else continue; 
      } else if (key[n]==POS) { if (vvo[i][j]<=0.) {fl=0; break;} else continue; 
      } else if (key[n]==NEG) { if (vvo[i][j]>=0.) {fl=0; break;} else continue; 
      } else if (key[n]==GTH) { if (vvo[i][j]<=arg[m]) {fl=0; break;} else continue; 
      } else if (key[n]==GTE) { if (vvo[i][j]< arg[m]) {fl=0; break;} else continue; 
      } else if (key[n]==LTH) { if (vvo[i][j]>=arg[m]) {fl=0; break;} else continue; 
      } else if (key[n]==LTE) { if (vvo[i][j]> arg[m]) {fl=0; break;} else continue; 
      } else if (key[n]==EQU) { if (vvo[i][j]!=arg[m]) {fl=0; break;} else continue; 
      } else if (key[n]==EQV) { if (vvo[i][j]!=vvo[i+1][j]) { fl=0; break;} else { 
          i++; continue; }
      } else if (key[n]==NEQ) { if (vvo[i][j]==arg[m]) {fl=0; break;} else continue; 
      } else if (key[n]==IBE) { if ((vvo[i][j]< arg[m])||(vvo[i][j]>=arg[m+1])) {
          fl=0; break; } else continue;  // IBE="[)" include-bracket-exclude
      } else if (key[n]==EBI) { if ((vvo[i][j]<=arg[m])||(vvo[i][j]> arg[m+1])) { 
          fl=0; break; } else continue;   // "(]" : exclude-bracket-include
      } else if (key[n]==IBI) { if ((vvo[i][j]< arg[m])||(vvo[i][j]> arg[m+1])) {
          fl=0; break; } else continue;   // "[]" : include-bracket-include
      } else if (key[n]==EBE) { if ((vvo[i][j]<=arg[m])||(vvo[i][j]>=arg[m+1])) {
          fl=0; break; } else continue;   // "()" : exclude-bracket-exclude
      } else {printf("vecst::slct ERR4 %g\n",key[n]); hoc_execerror("Unknown key",0);}
    }
    if (fl) ind[k++]=j; // all equal
  }
  vector_resize(vv, k);
  return k;
}
ENDVERBATIM

:* v.iwr()
: picks out numbers from multiple vectors using index ind
VERBATIM
static double iwr(void* vv) {
  int i, j, nx;
  double *x;
  FILE* f, *hoc_obj_file_arg();
  f = hoc_obj_file_arg(1);
  nx = vector_instance_px(vv, &x);
  if (nx>scrsz) { 
    if (scrsz>0) { free(scr); scr=(int *)NULL; }
    scrsz=nx+10000;
    scr=(int *)ecalloc(scrsz, sizeof(int));
  }
  for (i=0;i<nx;i++) scr[i]=(int)x[i]; /* copy into integer array */
  fwrite(&nx,sizeof(int),1,f);  // write out the size
  fwrite(scr,sizeof(int),nx,f);
  return nx;
}
ENDVERBATIM

:* v.ird()
: picks out numbers from multiple vectors using index ind
VERBATIM
static double ird(void* vv) {
  int i, j, nx, n;
  double *x;
  FILE* f, *hoc_obj_file_arg();
  f = hoc_obj_file_arg(1);
  nx = vector_instance_px(vv, &x);
  fread(&n,sizeof(int),1,f);  // size
  if (n>scrsz) { 
    if (scrsz>0) { free(scr); scr=(int *)NULL; }
    scrsz=n+10000;
    scr=(int *)ecalloc(scrsz, sizeof(int));
  }
  if (n!=nx) { 
    nx=vector_buffer_size(vv);
    if (n<=nx) {
      vector_resize(vv, n); nx=n; 
    } else {
      printf("%d > %d :: ",n,nx);
      hoc_execerror("Vector max capacity too small for ird ", 0);
    }
  }
  fread(scr,sizeof(int),n,f);
  for (i=0;i<nx;i++) x[i]=(double)scr[i];
  return n;
}
ENDVERBATIM

:* ind.insct(ind,vec1,vec2)
: picks out numbers from multiple vectors using index ind
VERBATIM
static double insct(void* vv) {
	int i, j, k, nx, nv1, nv2, maxsz;
	double *x, *v1, *v2;
	nx = vector_instance_px(vv, &x);
        maxsz=vector_buffer_size(vv);
        vector_resize(vv, maxsz);
	nv1 = vector_arg_px(1, &v1);
	nv2 = vector_arg_px(2, &v2);
        for (i=0,k=0;i<nv1;i++) for (j=0;j<nv2;j++) if (v1[i]==v2[j]) {
          if (k<maxsz) { x[k++]=v1[i]; } else {k++;}}  /* v1[i] found in both vectors */
        if (k>maxsz) { 
          printf("\tinsct WARNING: ran out of room: %d<%d\n",maxsz,k);
        } else { vector_resize(vv, k); }
	return (double)k;
}
ENDVERBATIM

:* vec.cull(src,key)
: remove numbers in vec that are found in the key
VERBATIM
static double cull(void* vv) {
	int i, j, k, nx, nv1, nv2, maxsz, flag;
	double *x, *v1, *v2;
	nx = vector_instance_px(vv, &x);
        maxsz=vector_buffer_size(vv);
        vector_resize(vv, maxsz);
	nv1 = vector_arg_px(1, &v1);
	nv2 = vector_arg_px(2, &v2);
        for (i=0,k=0;i<nv1;i++) {
          flag=1;
          for (j=0;j<nv2;j++) if (v1[i]==v2[j]) flag=0;
          if (flag) {if (k<maxsz) { x[k++]=v1[i]; } else { k++; }}
        }
        if (k>maxsz) { 
          printf("\tcull WARNING: ran out of room: %d<%d\n",maxsz,k);
        } else { vector_resize(vv, k); }
	return (double)k;
}
ENDVERBATIM

:* dest.redundout(src[,INDFLAG])
: remove numbers in vec that are found in the key
VERBATIM
static double redundout(void* vv) {
	int i, j, nx, nv1, maxsz, indflag;
	double *x, *v1, val;
        if (ifarg(2)) indflag=1; else indflag=0; 
	nx = vector_instance_px(vv, &x);
        maxsz=vector_buffer_size(vv);
        vector_resize(vv, maxsz);
	nv1 = vector_arg_px(1, &v1);
        val=v1[0]; x[0]=(indflag?0:val);
        for (j=1,i=1;i<nv1&&j<maxsz;i++) if (v1[i]!=val) { val=v1[i]; x[j++]=(indflag?i:val); }
        if (j==maxsz) { 
          printf("\tredundout WARNING: ran out of room: %d<needed\n",maxsz);
        } else { vector_resize(vv, j); }
	return (double)j;
}
ENDVERBATIM

:* vdest.cvlv(vsrc,vfilt)
: convolution
VERBATIM
static double cvlv(void* vv) {
  int i, j, k, nx, nsrc, nfilt;
  double *x, *src, *filt, sum, lpad, rpad;
  nx = vector_instance_px(vv, &x);
  nsrc = vector_arg_px(1, &src);
  nfilt = vector_arg_px(2, &filt);
  if (nx!=nsrc) { hoc_execerror("Vectors not same size: ", 0); }
  if (nfilt>nsrc) { hoc_execerror("Filter bigger than source ", 0); }
  for (i=0;i<nx;i++) {
    x[i]=0.0;
    for (j=0,k=i-(int)(nfilt/2);j<nfilt;j++,k++) {
      if (k>0 && k<nsrc-1) x[i]+=filt[j]*src[k];
    }
  }
}
ENDVERBATIM

:* vdest.intrp(flag)
: interpolate numbers replacing numbers given as flag
VERBATIM
static double intrp(void* vv) {
  int i, la, lb, nx;
  double *x, fl, a, b;
  nx = vector_instance_px(vv, &x);
  fl = *getarg(1);
  i=0; a=x[0]; la=0;
  if (a==fl) a=0;
  while (i<nx-1) {
    for (i=la+1;x[i]==fl && i<nx-1; i++) ; /* find the next one */
    b=x[i]; lb=i;
    for (i=la+1; i<lb; i++) x[i]= a + (b-a)/(lb-la)*(i-la);
    a=b; la=lb;
  }
}
ENDVERBATIM

:* tmp.nind(ind,vec1,vec2[,vec3,vec4,...])
: picks out numbers not in ind from multiple vectors
: ind must be sorted
VERBATIM
static double nind(void* vv) {
	int i, j, k, m, nx, ni, nv[10], num, c, last;
	double *x, *ind, *vvo[10];
	nx = vector_instance_px(vv, &x);
        for (i=0;ifarg(i);i++);
	if (i>9) hoc_execerror("ERR: nind can only handle 9 vectors", 0);
	num = i-2; /* number of vectors to be picked apart */
        for (i=0;i<num;i++) { 
          nv[i] = vector_arg_px(i+2, &vvo[i]);
          if (nx!=nv[i]) { printf("nind ERR %d %d %d\n",i,nx,nv[i]);
            hoc_execerror("Vectors must all be same size: ", 0); }
        }
        ni = vector_arg_px(1, &ind);
        c = nx-ni; /* the elems indexed are to be eliminated */
        if (ni>scrsz) { 
          if (scrsz>0) { free(scr); scr=(int *)NULL; }
          scrsz=ni+10000;
          scr=(int *)ecalloc(scrsz, sizeof(int));
        }
        for (i=0,last=-1;i<ni;i++) { 
          scr[i]=(int)ind[i]; /* copy into integer array */
          if (scr[i]<0 || scr[i]>=nx) hoc_execerror("nind(): Index out of bounds", 0);
          if (scr[i]<=last) hoc_execerror("nind(): indices should mono increase", 0);
          last=scr[i];
        }
        for (j=0;j<num;j++) { /* each output vec */
          for (i=0,last=-1,m=0;i<ni;i++) { /* the indices of ind */
            for (k=last+1;k<scr[i];k++) { x[m++]=vvo[j][k]; }
            last=scr[i];
          }
          for (k=last+1;k<nx;k++,m++) { x[m]=vvo[j][k]; }
          for (i=0;i<c;i++) vvo[j][i]=x[i];   
          vv=vector_arg(j+2); vector_resize(vv, c);
        }
	return c;
}
ENDVERBATIM

:* ind.keyind(key,vec1,vec2[,vec3,vec4,...])
: picks out indices of numbers in key from multiple vectors
VERBATIM
static double keyind(void* vv) {
  int i, j, k, ni, nk, nv[10], num;
  double *ind, *key, *vvo[10];
  ni = vector_instance_px(vv, &ind); // vv is ind
  for (i=0;ifarg(i);i++); i--; // drop back by one to get numarg()
  if (i>10) hoc_execerror("ERR: keyind can only handle 9 vectors", 0);
  num = i-1; /* number of vectors to be picked apart */
  for (i=0;i<num;i++) { 
    nv[i] = vector_arg_px(i+2, &vvo[i]);
    if (ni!=nv[i]) { printf("keyind ERR %d %d %d\n",i,ni,nv[i]);
    hoc_execerror("Non-key vectors must be same size: ", 0); }
  }
  nk = vector_arg_px(1, &key);
  if (nk!=num) { printf("keyind ERR2 %d %d\n",nk,num); 
    hoc_execerror("Key length must be number of vecs",0); }
  k=0;
  for (j=0;j<ni;j++) { // j steps through elements of vectors
    for (i=0;i<nk;i++) { // i steps through the key
      if (key[i]==ALL) continue; // OK - do nothing
      if (key[i]==NOZ)        { if (vvo[i][j]==0.) break; else continue; 
      } else if (key[i]==POS) { if (vvo[i][j]<=0.) break; else continue; 
      } else if (key[i]==NEG) { if (vvo[i][j]>=0.) break; else continue; 
      } else if (key[i]!=vvo[i][j]) break; // default
    }
    if (i==nk) ind[k++]=j; // all equal
  }
  vector_resize(vv, k);
  return k;
}
ENDVERBATIM

:* v1.flipbits(scratch,num) flips num bits
: uses scratch vector of same size as v1 to make sure doesn't flip same bit twice
VERBATIM
static double flipbits(void* vv) {	
	int i, nx, ny, flip, ii;
	double* x, *y;

	nx = vector_instance_px(vv, &x);
	ny = vector_arg_px(1, &y);
        flip = (int)*getarg(2);
	if (nx != ny) {
	  hoc_execerror("Scratch vector must be same size", 0);
	}
	for (i=0; i<nx; i++) { y[i]=x[i]; } /* copy */
	for (i=0; i < flip; i++) { /* flip these bits */
	  ii = (int) ((nx+1)*drand48());
	  if (x[ii]==y[ii]) { /* hasn't been touched */
	    x[ii]=((x[ii]==1.)?BVBASE:1.);
	  } else {
	    i--; /* do it again */
	  }
	}
	return flip;
}
ENDVERBATIM
 
:* v1.flipbalbits(scratch,num) flips num bits making sure to balance every 1
: flip with a 0 flip to preserve initial power
: uses scratch vector of same size as v1 to make sure doesn't flip same bit twice
VERBATIM
static double flipbalbits(void* vv) {	
	int i, nx, ny, flip, ii, next;
	double* x, *y;

	nx = vector_instance_px(vv, &x);
	ny = vector_arg_px(1, &y);
        flip = (int)*getarg(2);
	if (nx != ny) {
	  hoc_execerror("Scratch vector must be same size", 0);
	}
	for (i=0; i<nx; i++) { y[i]=x[i]; } /* copy */
        next = 1; /* start with 1 */
	for (i=0; i < flip;) { /* flip these bits */
	  ii = (int) ((nx+1)*drand48());
	  if (x[ii]==y[ii] && y[ii]==next) { /* hasn't been touched */
	    next=x[ii]=((x[ii]==1.)?BVBASE:1.);
            i++;
	  }
	}
	return flip;
}
ENDVERBATIM
 
:* v1.vpr() prints out neatly
VERBATIM
static double vpr(void* vv) {
  int i, nx;
  double* x;
  FILE* f;
  nx = vector_instance_px(vv, &x);
  if (ifarg(1)) { 
    f = hoc_obj_file_arg(1);
    for (i=0; i<nx; i++) {
      if (x[i]>BVBASE) { fprintf(f,"%d",1); 
      } else { fprintf(f,"%d",0); }
    }
    fprintf(f,"\n");
  } else {
    for (i=0; i<nx; i++) {
      if (x[i]>BVBASE) { printf("%d",1); 
      } else { printf("%d",0); }
    }
    printf("\n");
  }
  return 1.;
}
ENDVERBATIM
 
:* v1.thresh() threshold above and below thresh
VERBATIM
static double thresh(void* vv) {
  int i, nx, ny, cnt;
  double *x, *y, th;
  nx = vector_instance_px(vv, &x);
  cnt=0;
  if (hoc_is_object_arg(1)) { 
    ny = vector_arg_px(1, &y); th=0;
    if (nx!=ny) { hoc_execerror("Vector sizes don't match in thresh.", 0); }
    for (i=0; i<nx; i++) { if (x[i]>=y[i]) { x[i]= 1.; cnt++;} else { x[i]= BVBASE; } }
  } else { th = *getarg(1);
    for (i=0; i<nx; i++) { if (x[i] >= th) { x[i]= 1.; cnt++;} else { x[i]= BVBASE; } }
  }
  return cnt;
}
ENDVERBATIM
 
:* v1.triplet() return location of a triplet
VERBATIM
static double triplet(void* vv) {
  int i, nx;
  double *x, *y, a, b;
  nx = vector_instance_px(vv, &x);
  a = *getarg(1); b = *getarg(2);
  for (i=0; i<nx; i+=3) if (x[i]==a&&x[i+1]==b) break;
  if (i<nx) return (double)i; else return -1.;
}
ENDVERBATIM

:* state.onoff(volts,OBon,thresh,dur,refr)
:  looks at volts vector to decide if have rearch threshold thresh
:  OBon takes account of burst dur and refractory period
VERBATIM
static double onoff(void* vv) {
  int i, j, n, nv, non, nt, nd, nr, num;
  double *st, *vol, *obon, *thr, *dur, *refr;
  n = vector_instance_px(vv, &st);
  nv   = vector_arg_px(1, &vol);
  non  = vector_arg_px(2, &obon);
  nt   = vector_arg_px(3, &thr);
  nd   = vector_arg_px(4, &dur);
  nr   = vector_arg_px(5, &refr);
  if (n!=nv||n!=non||n!=nt||n!=nd||n!=nr) {
    hoc_execerror("v.onoff: vectors not all same size", 0); }
  for (i=0,num=0;i<n;i++) {
    obon[i]--;
    if (obon[i]>0.) { st[i]=1.; continue; } /* cell must fire */
    if (vol[i]>=thr[i] && obon[i]<= -refr[i]) { /* past refractory period */
        st[i]=1.; obon[i]=dur[i]; num++;
    } else { st[i]= BVBASE; }
  }
  return num;
}
ENDVERBATIM
 
:* vo.bpeval(outp,del)
:  service routine for back-prop: vo=outp*(1-outp)*del
VERBATIM
static double bpeval(void* vv) {
  int i, n, no, nd, flag=0;
  double add,div;
  double *vo, *outp, *del;
  n = vector_instance_px(vv, &vo);
  no   = vector_arg_px(1, &outp);
  nd   = vector_arg_px(2, &del);
  if (ifarg(3) && ifarg(4)) { add= *getarg(3); div= *getarg(4); flag=1;}
  if (n!=no||n!=nd) hoc_execerror("v.bpeval: vectors not all same size", 0);
  if (flag) {
    for (i=0;i<n;i++) vo[i]=((outp[i]+add)/div)*(1.-1.*((outp[i]+add)/div))*del[i];
  } else {
    for (i=0;i<n;i++) vo[i]=outp[i]*(1.-1.*outp[i])*del[i];
  }
}
ENDVERBATIM
 
:* v1.sedit() string edit
VERBATIM
static double sedit(void* vv) {
  int i, n, ni, f=0;
  double *x, *ind, th, val;
  Symbol* s; char *op;
  op = gargstr(1);
  n = vector_instance_px(vv, &x);
  sprintf(op,"hello world");
  return (double)n;
}
ENDVERBATIM

:* v1.w() a .where that sets elements in source vector
VERBATIM
static double w(void* vv) {
  int i, n, ni, f=0;
  double *x, *ind, th, val;
  Symbol* s; char *op;
  if (! ifarg(1)) { 
    printf("v1.w('op',thresh[,val,v2])\n"); 
    printf("  a .where that sets elements in v1 to val (default 0), if v2 => only look at these elements\n");
    printf("  'op'=='function name' is a .apply targeted by v2 called as func(x[i],thresh,val)\n");
    return -1.;
  }
  op = gargstr(1);
  n = vector_instance_px(vv, &x);
  th = *getarg(2);
  if (ifarg(3)) { val = *getarg(3); } else { val = 0.0; }
  if (ifarg(4)) {ni = vector_arg_px(4, &ind); f=1;} // just look at the spots indexed
  if (!strcmp(op,"==")) { 
    if (f==1) {for (i=0; i<ni;i++) {if (x[(int)ind[i]]==th) { x[(int)ind[i]] = val;}}
    } else {for (i=0; i<n; i++) {if (x[i]==th) { x[i] = val;}}}
  } else if (!strcmp(op,"!=")) {
    if (f==1) {for (i=0; i<ni;i++) {if (x[(int)ind[i]]!=th) { x[(int)ind[i]] = val;}}
    } else {for (i=0; i<n; i++) {if (x[i]!=th) { x[i] = val;}}}
  } else if (!strcmp(op,">")) {
    if (f==1) {for (i=0; i<ni;i++) {if (x[(int)ind[i]]>th) { x[(int)ind[i]] = val;}}
    } else {for (i=0; i<n; i++) {if (x[i]>th) { x[i] = val;}}}
  } else if (!strcmp(op,"<")) {
    if (f==1) {for (i=0; i<ni;i++) {if (x[(int)ind[i]]<th) { x[(int)ind[i]] = val;}}
    } else {for (i=0; i<n; i++) {if (x[i]<th) { x[i] = val;}}}
  } else if (!strcmp(op,">=")) {
    if (f==1) {for (i=0; i<ni;i++) {if (x[(int)ind[i]]>=th) { x[(int)ind[i]] = val;}}
    } else {for (i=0; i<n; i++) {if (x[i]>=th) { x[i] = val;}}}
  } else if (!strcmp(op,"<=")) {
    if (f==1) {for (i=0; i<ni;i++) {if (x[(int)ind[i]]<=th) { x[(int)ind[i]] = val;}}
    } else {for (i=0; i<n; i++) {if (x[i]<=th) { x[i] = val;}}}
  } else if ((s=hoc_lookup(op))) { // same as .apply but only does indexed ones
    if (f==1) {for (i=0; i<ni;i++) {
        hoc_pushx(x[(int)ind[i]]); hoc_pushx(th); hoc_pushx(val);
      x[(int)ind[i]]=hoc_call_func(s, 3);}
    } else {for (i=0; i<n;i++) {hoc_pushx(x[i]); hoc_pushx(th); hoc_pushx(val);
        x[i]=hoc_call_func(s, 3);}}
  }
  return (double)i;
}
ENDVERBATIM
  
:* dest.xing(src,tvec,thresh) 
: a .where that looks for threshold crossings and then doesn't pick another till
: comes back down again; places values from tvec in dest; interpolates
VERBATIM
static double xing(void* vv) {
  int i, j, nsrc, ndest, ntvec, f, maxsz;
  double *src, *dest, *tvec, th;
  ndest = vector_instance_px(vv, &dest);
  nsrc = vector_arg_px(1, &src);
  ntvec = vector_arg_px(2, &tvec);
  th = *getarg(3);
  maxsz=vector_buffer_size(vv);
  vector_resize(vv, maxsz);
  if (nsrc!=ntvec) hoc_execerror("v.xing: vectors not all same size", 0);
  for (i=0,f=0,j=0; i<nsrc; i++) {
    if (src[i]>th) { /* ? passing thresh */
      if (f==0) { 
        if (j>=maxsz) {
          printf("(%d) :: ",maxsz);
          hoc_execerror("Dest vec too small in xing ", 0);
        }
        if (i>0) { /* don't record if first value is above thresh */
          dest[j++] = tvec[i-1] + (tvec[i]-tvec[i-1])*(th-src[i-1])/(src[i]-src[i-1]);
        }
        f=1; 
      }
    } else {       /* below thresh */
      if (f==1) { f=0; } /* just passed going down */
    }
  }
  vector_resize(vv, j);
  return (double)i;
}
ENDVERBATIM

:* v1.xzero() looks for zero crossings 
VERBATIM
static double xzero(void* vv) {
  int i, n, nv, up;
  double *x, *vc, th, cnt=0.;
  n = vector_instance_px(vv, &x);
  nv = vector_arg_px(1, &vc);
  if (ifarg(2)) { th = *getarg(2); } else { th=0.0;}
  if (vc[0]<th) up=0; else up=1;  /* F or T */
  if (nv!=n) hoc_execerror("Vector size doesn't match.", 0);
  for (i=0; i<nv; i++) {
    x[i]=0.;
    if (up) { /* look for passing down */
      if (vc[i]<th) { x[i]=-1; up=0; cnt++; }
    } else if (vc[i]>th) { up=x[i]=1; cnt++; }
  }
  return cnt;
}
ENDVERBATIM
  
:* v1.sw(FROM,TO) switchs all FROMs to TO
VERBATIM
static double sw(void* vv) {
  int i, n;
  double *x, fr, to;
  n = vector_instance_px(vv, &x);
  fr = *getarg(1);
  to = *getarg(2);
  for (i=0; i<n; i++) {
    if (x[i]==fr) { x[i] = to;}
  }
  return n;
}
ENDVERBATIM
 
:* v.v2d(&x) copies from vector to double area -- a seg error waiting to happen
VERBATIM
static double v2d(void* vv) {
  int i, n, num;
  double *x, *to;
  n = vector_instance_px(vv, &x);
  to = hoc_pgetarg(1);
  if (ifarg(2)) { num = *getarg(2); } else { num=-1;}
  if (num>-1 && num!=n) { hoc_execerror("Vector size doesn't match.", 0); }
  for (i=0; i<n; i++) {to[i] = x[i];}
  return n;
}
ENDVERBATIM
 
:* v.d2v(&x) copies from double area to vector -- a seg error waiting to happen
VERBATIM
static double d2v(void* vv) {
  int i, n, num;
  double *x, *fr;
  n = vector_instance_px(vv, &x);
  fr = hoc_pgetarg(1);
  if (ifarg(2)) { num = *getarg(2); } else { num=-1;}
  if (num>-1 && num!=n) { hoc_execerror("Vector size doesn't match.", 0); }
  for (i=0; i<n; i++) {x[i] = fr[i];}
  return n;
}
ENDVERBATIM
 
:* v1.ismono([arg]) asks whether is monotonically increasing, with arg==-1 - decreasing
:  with arg==0 - all same
VERBATIM
static double ismono(void* vv) {
  int i, n, flag;
  double *x,last;
  n = vector_instance_px(vv, &x);
  if (ifarg(1)) { flag = (int)*getarg(1); } else { flag = 1; }
  last=x[0]; 
  if (flag==1) {
    for (i=1; i<n && x[i]>=last; i++) last=x[i];
  } else if (flag==-1) {
    for (i=1; i<n && x[i]<=last; i++) last=x[i];
  } else if (flag==0) {
    for (i=1; i<n && x[i]==last; i++) ;
  }
  if (i==n) return 1.; else return 0.;
}
ENDVERBATIM
 
:* v1.count(num) returns number of instances of num
VERBATIM
static double count(void* vv) {
  int i, n, cnt=0;
  double *x,num;
  n = vector_instance_px(vv, &x);
  num = *getarg(1);
  for (i=0; i<n; i++) if (x[i]==num) cnt++;
  return cnt;
}
ENDVERBATIM

:* fac (n) 
: from numerical recipes p.214
FUNCTION fac (n) {
VERBATIM {
    static int ntop=4;
    static double a[101]={1.,1.,2.,6.,24.};
    static double cof[6]={76.18009173,-86.50532033,24.01409822,
      -1.231739516,0.120858003e-2,-0.536382e-5};
    int j,n;
    n = (int)_ln;
    if (n<0) { hoc_execerror("No negative numbers ", 0); }
    if (n>100) { /* gamma function */
      double x,tmp,ser;
      x = _ln;
      tmp=x+5.5;
      tmp -= (x+0.5)*log(tmp);
      ser=1.0;
      for (j=0;j<=5;j++) {
        x += 1.0;
        ser += cof[j]/x;
      }
      return exp(-tmp+log(2.50662827465*ser));
    } else {
      while (ntop<n) {
        j=ntop++;
        a[ntop]=a[j]*ntop;
      }
    return a[n];
    }
}
ENDVERBATIM
}
 
:* logfac (n)
: from numerical recipes p.214
FUNCTION logfac (n) {
VERBATIM {
    static int ntop=4;
    static double a[101]={1.,1.,2.,6.,24.};
    static double cof[6]={76.18009173,-86.50532033,24.01409822,
      -1.231739516,0.120858003e-2,-0.536382e-5};
    int j,n;
    n = (int)_ln;
    if (n<0) { hoc_execerror("No negative numbers ", 0); }
    if (n>100) { /* gamma function */
      double x,tmp,ser;
      x = _ln;
      tmp=x+5.5;
      tmp -= (x+0.5)*log(tmp);
      ser=1.0;
      for (j=0;j<=5;j++) {
        x += 1.0;
        ser += cof[j]/x;
      }
      return (-tmp+log(2.50662827465*ser));
    } else {
      while (ntop<n) {
        j=ntop++;
        a[ntop]=a[j]*ntop;
      }
    return log(a[n]);
    }
}
ENDVERBATIM
}

PROCEDURE Expo (x) {
  TABLE RES FROM -20 TO 0 WITH 5000
  RES = exp(x)
}

FUNCTION AAA (x) {
  Expo(x)
  AAA = RES
}

:* dest.smgs(src,low,high,step,var)
:  rewrite on v.sumgauss() in nrn5.3::ivoc/ivocvect.cpp:1078
VERBATIM
static double smgs(void* vv) {	
  int i, j, nx, xv, nsum, points, maxsz;
  double *x, *sum;
  double  low , high , step , var , svar , scale , arg;

  nsum = vector_instance_px(vv, &sum);
  nx = vector_arg_px(1,&x);
  low = *getarg(2);
  high = *getarg(3);
  step = *getarg(4);
  var = *getarg(5);

  points = (int)((high-low)/step+hoc_epsilon);
  if (nsum!=points) { 
    maxsz=vector_buffer_size(vv);
    if (points<=maxsz) {
      vector_resize(vv, nsum); nsum=points; 
    } else {
      printf("%d > %d :: ",points,maxsz);
      hoc_execerror("Vector max capacity too small in smgs ", 0);
    }
  }

  svar = -2.*var*var/step/step;
  scale = 1./sqrt(2.*M_PI)/var;

  for (j=0; j<points;j++) sum[j] = 0.;
  for (i=0;i<nx;i++) {
    xv = (int)((x[i]-low)/step + 0.5);
    for (j=xv; j<points && (arg=(j-xv)*(j-xv)/svar)>-20;j++) {
      Expo(arg);
      sum[j] += RES;
    }
    for (j=xv-1; j>=0 && (arg=(j-xv)*(j-xv)/svar)>-20;j--) {
      Expo(arg);
      sum[j] += RES;
    }
  }
  for (j=0; j<points;j++) sum[j] *= scale;
  return svar;
}
ENDVERBATIM

:* PROCEDURE install_vecst()
PROCEDURE install_vecst() {
  VECST_INSTALLED=1
VERBATIM
  install_vector_method("slope", slope);
  install_vector_method("vslope", vslope);
  install_vector_method("stats", stats);
  install_vector_method("vstats", vstats);
  install_vector_method("indset", indset);
  install_vector_method("randwd", randwd);
  install_vector_method("hamming", hamming);
  install_vector_method("flipbits", flipbits);
  install_vector_method("flipbalbits", flipbalbits);
  install_vector_method("vpr", vpr);
  install_vector_method("thresh", thresh);
  install_vector_method("triplet", triplet);
  install_vector_method("onoff", onoff);
  install_vector_method("bpeval", bpeval);
  install_vector_method("w", w);
  install_vector_method("sedit", sedit);
  install_vector_method("xing", xing);
  install_vector_method("cvlv", cvlv);
  install_vector_method("intrp", intrp);
  install_vector_method("xzero", xzero);
  install_vector_method("sw", sw);
  install_vector_method("ismono", ismono);
  install_vector_method("count", count);
  install_vector_method("fewind", fewind);
  install_vector_method("findx", findx);
  install_vector_method("sindx", sindx);
  install_vector_method("sindv", sindv);
  install_vector_method("nind", nind);
  install_vector_method("keyind", keyind);
  install_vector_method("slct", slct);
  install_vector_method("insct", insct);
  install_vector_method("cull", cull);
  install_vector_method("redundout", redundout);
  install_vector_method("d2v", d2v);
  install_vector_method("v2d", v2d);
  install_vector_method("iwr", iwr);
  install_vector_method("ird", ird);
  install_vector_method("smgs", smgs);
ENDVERBATIM
}

: unable to get the drand here to recognize the same fseed used in rand
PROCEDURE vseed(seed) {
  VERBATIM
  srand48((unsigned)_lseed);
  ENDVERBATIM
}

Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ (1998) Computational models of thalamocortical augmenting responses. J Neurosci 18:6444-65[PubMed]

References and models cited by this paper

References and models that cite this paper

Avoli M (1986) Inhibitory potentials in neurons of the deep layers of the in vitro neocortical slice. Brain Res 370:165-70

Bal T, McCormick DA (1996) What stops synchronized thalamocortical oscillations? Neuron 17:297-308 [PubMed]

Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ (1998) Cellular and network models for intrathalamic augmenting responses during 10-Hz stimulation. J Neurophysiol 79:2730-48 [Journal] [PubMed]

Castro-Alamancos MA, Connors BW (1996) Spatiotemporal properties of short-term plasticity sensorimotor thalamocortical pathways of the rat. J Neurosci 16:2767-79 [PubMed]

Castro-Alamancos MA, Connors BW (1996) Cellular mechanisms of the augmenting response: short-term plasticity in a thalamocortical pathway. J Neurosci 16:7742-56 [PubMed]

Connors BW, Malenka RC, Silva LR (1988) Two inhibitory postsynaptic potentials, and GABAA and GABAB receptor-mediated responses in neocortex of rat and cat. J Physiol 406:443-68 [PubMed]

Contreras D, Curro Dossi R, Steriade M (1993) Electrophysiological properties of cat reticular thalamic neurones in vivo. J Physiol 470:273-94 [PubMed]

Contreras D, Destexhe A, Sejnowski TJ, Steriade M (1997) Spatiotemporal patterns of spindle oscillations in cortex and thalamus. J Neurosci 17:1179-96 [PubMed]

Contreras D, Steriade M (1995) Cellular basis of EEG slow rhythms: a study of dynamic corticothalamic relationships. J Neurosci 15:604-22 [PubMed]

Contreras D, Steriade M (1996) Spindle oscillation in cats: the role of corticothalamic feedback in a thalamically generated rhythm. J Physiol 490 ( Pt 1):159-79 [PubMed]

Contreras D, Timofeev I, Steriade M (1996) Mechanisms of long-lasting hyperpolarizations underlying slow sleep oscillations in cat corticothalamic networks. J Physiol 494 ( Pt 1):251-64 [PubMed]

Crunelli V, Haby M, Jassik-Gerschenfeld D, Leresche N, Pirchio M (1988) Cl- - and K+-dependent inhibitory postsynaptic potentials evoked by interneurones of the rat lateral geniculate nucleus. J Physiol 399:153-76 [PubMed]

Destexhe A, Bal T, McCormick DA, Sejnowski TJ (1996) Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J Neurophysiol 76:2049-70 [Journal] [PubMed]

   Thalamocortical and Thalamic Reticular Network (Destexhe et al 1996) [Model]

Destexhe A, Contreras D, Sejnowski TJ, Steriade M (1994) A model of spindle rhythmicity in the isolated thalamic reticular nucleus. J Neurophysiol 72:803-18 [Journal] [PubMed]

   Thalamic Reticular Network (Destexhe et al 1994) [Model]

Destexhe A, Contreras D, Steriade M, Sejnowski TJ, Huguenard JR (1996) In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. J Neurosci 16:169-85 [Journal] [PubMed]

   [1 reconstructed morphology on NeuroMorpho.Org]
   Thalamic reticular neurons: the role of Ca currents (Destexhe et al 1996) [Model]

Destexhe A, Mainen ZF, Sejnowski TJ (1994) Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism. J Comput Neurosci 1:195-230 [Journal] [PubMed]

   Application of a common kinetic formalism for synaptic models (Destexhe et al 1994) [Model]
   Kinetic synaptic models applicable to building networks (Destexhe et al 1998) [Model]

Durmuller N, Timofeev I, Grenier F, Steriade M (1998) Mechanisms of intracortical augmenting responses Soc Neurosci Abstr

Dutar P, Nicoll RA (1988) A physiological role for GABAB receptors in the central nervous system. Nature 332:156-8 [PubMed]

Enright WH, Higham DJ, Owren B, Sharp PW (1995) A survey of the explicit Runge-Kutta method. Available from ftp://ftp.cs.toronto.edu/pub/reports/na/cs-94-291.ps.Z

Ferster D, Lindstrom S (1985) Augmenting responses evoked in area 17 of the cat by intracortical axon collaterals of cortico-geniculate cells. J Physiol 367:217-32 [PubMed]

Hernandez-Cruz A, Pape HC (1989) Identification of two calcium currents in acutely dissociated neurons from the rat lateral geniculate nucleus. J Neurophysiol 61:1270-83 [Journal] [PubMed]

Hirsch JC, Burnod Y (1987) A synaptically evoked late hyperpolarization in the rat dorsolateral geniculate neurons in vitro. Neuroscience 23:457-68 [PubMed]

Hodgkin AL, Huxley AF (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 117:500-44 [Journal] [PubMed]

   Squid axon (Hodgkin, Huxley 1952) (LabAXON) [Model]
   Squid axon (Hodgkin, Huxley 1952) (NEURON) [Model]
   Squid axon (Hodgkin, Huxley 1952) (SNNAP) [Model]
   Squid axon (Hodgkin, Huxley 1952) used in (Chen et al 2010) (R language) [Model]
   Squid axon (Hodgkin, Huxley 1952) (SBML, XPP, other) [Model]

Houweling AR, Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ (1998) Computational models of intracortical augmenting responses Soc Neurosci Abstr

Huguenard JR, Coulter DA, Prince DA (1991) A fast transient potassium current in thalamic relay neurons: kinetics of activation and inactivation. J Neurophysiol 66:1304-15 [Journal] [PubMed]

Huguenard JR, McCormick DA (1992) Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J Neurophysiol 68:1373-83 [Journal] [PubMed]

Huguenard JR, Prince DA (1992) A novel T-type current underlies prolonged Ca(2+)-dependent burst firing in GABAergic neurons of rat thalamic reticular nucleus. J Neurosci 12:3804-17 [PubMed]

Huguenard JR, Prince DA (1994) Intrathalamic rhythmicity studied in vitro: nominal T-current modulation causes robust antioscillatory effects. J Neurosci 14:5485-502 [PubMed]

Jahnsen H, Llinas R (1984) Ionic basis for the electro-responsiveness and oscillatory properties of guinea-pig thalamic neurones in vitro. J Physiol 349:227-47 [PubMed]

Jahnsen H, Llinas R (1984) Electrophysiological properties of guinea-pig thalamic neurones: an in vitro study. J Physiol 349:205-26 [PubMed]

Jones EG (1985) The Thalamus

Kammermeier PJ, Jones SW (1997) High-voltage-activated calcium currents in neurons acutely isolated from the ventrobasal nucleus of the rat thalamus. J Neurophysiol 77:465-75 [Journal] [PubMed]

Kandel A, Buzsaki G (1997) Cellular-synaptic generation of sleep spindles, spike-and-wave discharges, and evoked thalamocortical responses in the neocortex of the rat. J Neurosci 17:6783-97 [PubMed]

Leresche N, Lightowler S, Soltesz I, Jassik-Gerschenfeld D, Crunelli V (1991) Low-frequency oscillatory activities intrinsic to rat and cat thalamocortical cells. J Physiol 441:155-74 [PubMed]

Luthi A, Mccormick DA (1997) Both electrophysiological and biochemical oscillations determine spindle wave rhythmicity Soc Neurosci Abstr 23:1820

Mainen ZF, Sejnowski TJ (1996) Influence of dendritic structure on firing pattern in model neocortical neurons. Nature 382:363-6 [Journal] [PubMed]

   [2 reconstructed morphologies on NeuroMorpho.Org]
   Pyramidal Neuron Deep, Superficial; Aspiny, Stellate (Mainen and Sejnowski 1996) [Model]

McCormick DA, Huguenard JR (1992) A model of the electrophysiological properties of thalamocortical relay neurons. J Neurophysiol 68:1384-400 [Journal] [PubMed]

McCormick DA, Pape HC (1990) Properties of a hyperpolarization-activated cation current and its role in rhythmic oscillation in thalamic relay neurones. J Physiol 431:291-318 [PubMed]

   Thalamic Relay Neuron: I-h (McCormick, Pape 1990) [Model]

Metherate R, Ashe JH (1994) Facilitation of an NMDA receptor-mediated EPSP by paired-pulse stimulation in rat neocortex via depression of GABAergic IPSPs. J Physiol 481 ( Pt 2):331-48 [PubMed]

Morin D, Steriade M (1981) Development from primary to augmenting responses in the somatosensory system. Brain Res 205:49-66 [PubMed]

Morison RS, Dempsey EW (1943) Mechanisms of thalamocortical augmentation and repetition Am J Physiol 138:297-308

Pare D, Dossi RC, Steriade M (1991) Three types of inhibitory postsynaptic potentials generated by interneurons in the anterior thalamic complex of cat. J Neurophysiol 66:1190-204 [Journal] [PubMed]

Pedroarena C, Llinas R (1997) Dendritic calcium conductances generate high-frequency oscillation in thalamocortical neurons. Proc Natl Acad Sci U S A 94:724-8

Soltesz I, Lightowler S, Leresche N, Jassik-Gerschenfeld D, Pollard CE, Crunelli V (1991) Two inward currents and the transformation of low-frequency oscillations of rat and cat thalamocortical cells. J Physiol 441:175-97 [PubMed]

Steriade M (1984) The excitatory-inhibitory response sequence of thalamic and neocortical cells: state related changes and regulatory systems Dynamic Aspects of Neocortical Function, Edelman GM:Gall WE:Cowan WM, ed. pp.107

Steriade M, Contreras D, Curro Dossi R, Nunez A (1993) 1 Hz) oscillation in reticular thalamic and thalamocortical neurons: scenario of sleep rhythm generation in interacting thalamic and neocortical networks. J Neurosci 13:3284-99

Steriade M, Domich L, Oakson G (1986) Reticularis thalami neurons revisited: activity changes during shifts in states of vigilance. J Neurosci 6:68-81 [PubMed]

Steriade M, Dossi RC, Nunez A (1991) Network modulation of a slow intrinsic oscillation of cat thalamocortical neurons implicated in sleep delta waves: cortically induced synchronization and brainstem cholinergic suppression. J Neurosci 11:3200-17 [PubMed]

Steriade M, Nunez A, Amzica F (1993) 1 Hz) oscillation of neocortical neurons in vivo: depolarizing and hyperpolarizing components. J Neurosci 13:3252-65 [PubMed]

Steriade M, Nunez A, Amzica F (1993) 1 Hz) neocortical oscillation and other sleep rhythms of the electroencephalogram. J Neurosci 13:3266-83 [PubMed]

Steriade M, Timofeev I (1997) Short-term plasticity during intrathalamic augmenting responses in decorticated cats. J Neurosci 17:3778-95 [PubMed]

Steriade M, Timofeev I, Grenier F, Durmuller N (1998) Role of thalamic and cortical neurons in augmenting responses and self-sustained activity: dual intracellular recordings in vivo. J Neurosci 18:6425-43 [PubMed]

Steriade M, Wyzinski P, Apostol V (1972) Corticofugal projections governing rhythmic thalamic activity Corticothalamic Projections And Sensorimotor Activities, Frigyesi TL:Rinvik E:Yahr MD, ed. pp.221

Timofeev I, Contreras D, Steriade M (1996) Synaptic responsiveness of cortical and thalamic neurones during various phases of slow sleep oscillation in cat. J Physiol 494 ( Pt 1):265-78 [PubMed]

Timofeev I, Steriade M (1996) Low-frequency rhythms in the thalamus of intact-cortex and decorticated cats. J Neurophysiol 76:4152-68 [Journal] [PubMed]

Timofeev I, Steriade M (1997) Fast (mainly 30-100 Hz) oscillations in the cat cerebellothalamic pathway and their synchronization with cortical potentials. J Physiol 504 ( Pt 1):153-68 [PubMed]

Timofeev I, Steriade M (1998) Cellular mechanisms underlying intrathalamic augmenting responses of reticular and relay neurons. J Neurophysiol 79:2716-29 [Journal] [PubMed]

Traub RD, Miles R (1991) Neuronal Networks Of The Hippocampus

Ulrich D, Huguenard JR (1997) Nucleus-specific chloride homeostasis in rat thalamus. J Neurosci 17:2348-54 [PubMed]

Zhou Q, Godwin DW, O'Malley DM, Adams PR (1997) Visualization of calcium influx through channels that shape the burst and tonic firing modes of thalamic relay cells. J Neurophysiol 77:2816-25 [Journal] [PubMed]

Chen XD, Shu SF, Kennedy DP, Willcox SC, Bayliss DA (2009) Subunit-specific effects of isoflurane on neuronal I-h in HCN1 knockout mice. J. Neurophysiol. 101:129-149 [Journal]

Houweling AR, Bazhenov M, Timofeev I, Grenier F, Steriade M, Sejnowski TJ (2002) Frequency-selective augmenting responses by short-term synaptic depression in cat neocortex. J Physiol 542:599-617 [PubMed]

Houweling AR, Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ (1999) Cortical and thalamic components of augmenting responses: A modeling study Neurocomputing 26-27:735-742

Karameh FN, Dahleh MA, Brown EN, Massaquoi SG (2006) Modeling the contribution of lamina 5 neuronal and network dynamics to low frequency EEG phenomena. Biol Cybern 95:289-310 [PubMed]

Li X, Morita K, Robinson HP, Small M (2013) Control of layer 5 pyramidal cell spiking by oscillatory inhibition in the distal apical dendrites: a computational modelling study. J Neurophysiol 109(11):2739-2756 [Journal] [PubMed]

   L5 pyr. cell spiking control by oscillatory inhibition in distal apical dendrites (Li et al 2013) [Model]

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, Hines ML (2005) Independent variable time-step integration of individual neurons for network simulations. Neural Comput 17:903-21 [Journal] [PubMed]

   Local variable time step method (Lytton, Hines 2005) [Model]

Rulkov NF, Timofeev I, Bazhenov M (2004) Oscillations in large-scale cortical networks: map-based model. J Comput Neurosci 17:203-23 [Journal] [PubMed]

   Large cortex model with map-based neurons (Rulkov et al 2004) [Model]

Sejnowski TJ, Destexhe A (2000) Why do we sleep? Brain Res 886:208-223 [PubMed]

Skorheim SW, Razak K, Bazhenov M (2014) Network models of frequency modulated sweep detection PLoS ONE 9(12):e115196 [Journal] [PubMed]

   Network models of frequency modulated sweep detection (Skorheim et al. 2014) [Model]

Sohal VS, Huntsman MM, Huguenard JR (2000) Reciprocal inhibitory connections regulate the spatiotemporal properties of intrathalamic oscillations. J Neurosci 20:1735-45 [PubMed]

(71 refs)