Auditory nerve model for predicting performance limits (Heinz et al 2001)

 Download zip file 
Help downloading and running models
Accession:36834
A computational auditory nerve (AN) model was developed for use in modeling psychophysical experiments with normal and impaired human listeners. In this phenomenological model, many physiologically vulnerable response properties associated with the cochlear amplifier are represented by a single nonlinear control mechanism, see paper for details. Several model versions are described that can be used to evaluate the relative effects of these nonlinear properties.
Reference:
1 . Heinz MG, Zhang X, Bruce IC, Carney LH (2001) Auditory nerve model for predicting performance limits of normal and impaired listeners. Acoustics Research Letters Online 2(3):91-96
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s): Auditory nerve;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program; MATLAB;
Model Concept(s): Activity Patterns; Temporal Pattern Generation; Audition;
Implementer(s): Zhang, Xuedong ;
// To do : 1. change the option that about # of model
//         2. add the option that can generate the synapse output / spike rate
//         3. check the parameter for CAT/HUMAN, and maybe more
//         4. Species can be specified as others as the same in the paper

 /* Model numbers as used in ARLO Heinz et al., Fig. 4 */
//  model == 1 bmmodel = FeedForward_NL;
//  model == 2 bmmodel = FeedBack_NL;
//  model == 3 bmmodel = Sharp_Linear;
//  model == 4 bmmodel = Broad_Linear;
//  model == 5 bmmodel = Broad_Linear_High;

#include "cmpa.h"
#include "hc.h"
#include "synapse.h"
#include "complex.h"
#include "filters.h"
#include <stdlib.h>
double Get_tau(int species, double cf, int order, double* taumax, double* taumin, double* taurange);
double erbGM(double);
double cochlea_f2x(int species,double f);
double cochlea_x2f(int species,double x);

void runAN2(TAuditoryNerve *p,const double *in,double *out,const int length);
double runBasilarMembrane(TBasilarMembrane *bm, double x);
void run2BasilarMembrane(TBasilarMembrane *bm, const double *in, double *out, const int length);
/*
 * different species have different parameters !!!!
 * 1. Get_Tau(...)
 * 2. Synapse setup in this function
 * 3. IHC low pass filter cutoff freq
 * 
 * */
void initAuditoryNerve(TAuditoryNerve *p,int model, int species, double tdres, double cf, double spont)
{ /* default */
  Tsynapse *syn;

  double Kcf,temp;
  double p1,p2,pst,psl;

  p->run = runAN;
  p->run2 = runAN2;

  syn = &(p->syn);

  p->model = model;
  p->species = species;
  p->cf = cf;
  p->spont = spont;
  p->tdres = tdres;

  p->run = runAN;
  p->run2 = runAN2;
  /*
   * Init the basilar membrane
   * */
  initBasilarMembrane(&(p->bm),model,species,tdres,cf);

  /* Now add support for different species */

  /* This line outputs the appropriate sout rate based on whether or not spikes are being generated (with refractoriness) */
  if(p->ifspike==1) syn->Ass = 350;
  else syn->Ass = 130; /* borrow this for Mike G. Heinz */
  switch(species)
  { 
	  case 1 : /* Cat - low CFs only (Carney '93) */
	  case 9 : /* Cat - all CFs "Universal"  (Zhang et al 2001) */
		  /*
		   * Init the synapse
		   * */
		  syn->tdres = tdres;
		  syn->cf = cf;
		  syn->spont = spont;
		  syn->Pimax = 0.6;
		  syn->PTS = 1+9*spont/(9+spont);
		  syn->Ar_over_Ast = 6.0;
		  syn->tauR = 0.002;
		  syn->tauST = 0.060;
		  initSynapse(syn);  /*This function calcuate all other parameters as described in the appendix*/

		  /*
		   * Init the ihcppi rectify function and parameters
		   * !!!!Use the result from initSynapse(syn)
		   *
		   * */
 	  	  Kcf = 2.0+3.0*log10(cf/1000);
		  if(Kcf<1.5) Kcf = 1.5;
		  syn->Vsat = syn->Pimax*Kcf*20.0*(1+spont)/(5+spont);
		  /* Important !!! The original code seems use the following equation!!!
		     This is not shown in the paper but Prest is very small compare to the other part < 0.01*Vsat
		     and the results is the same
		   */
		  /* syn->Vsat = syn->Pimax*Kcf*20.0*(1+spont)/(5+spont)+syn->Prest; */
		  
		  break;
	  case 0 :
	  default: /* Human */
		  /*
		   * Init the synapse
		   * */
		  syn->tdres = tdres;
		  syn->cf = cf;
		  syn->spont = spont;
		  syn->Pimax = 0.6;
		  syn->PTS = 1+9*spont/(9+spont);   /* Peak to Steady State Ratio, characteristic of PSTH */
		  syn->Ar_over_Ast = 6.0;
		  syn->tauR = 0.002;
		  syn->tauST = 0.060;
		  initSynapse(syn);  /*This function calcuate all other parameters as described in the appendix*/

		  /*
		   * Init the ihcppi rectify function and parameters
		   * !!!!Use the result from initSynapse(syn)
		   *
		   * */
            	  Kcf = 2.0+1.3*log10(cf/1000); /*/ From M. G. Heinz */
		  if(Kcf<1.5) Kcf = 1.5;
		  syn->Vsat = syn->Pimax*Kcf*60.0*(1+spont)/(6+spont);
		  /* Important !!! The original code seems use the following equation!!!
		     This is not shown in the paper but Prest is very small compare to the other part < 0.01*Vsat
		     and the results is the same, so we use the upper equation
		   */
		  /* syn->Vsat = syn->Pimax*Kcf*60.0*(1+spont)/(6+spont)+syn->Prest; */
		  
		  break;
  };
  temp = log(2)*syn->Vsat/syn->Prest;
  if(temp<400) pst = log(exp(temp)-1);
  else pst = temp;
  psl = syn->Prest*pst/log(2);
  p2 = pst;
  p1 = psl/pst;

#ifdef DEBUG
  printf("\ninitSyanpse : p1=%f; p2=%f",p1,p2);
#endif
  
  p->ihcppi.psl = psl;
  p->ihcppi.pst = pst;
  p->ihcppi.p1 = p1;
  p->ihcppi.p2 = p2;
  p->ihcppi.run = runIHCPPI;	/*These are functions used to get ihcppi output */
  p->ihcppi.run2 = runIHCPPI2;  /*Defined in hc.c, hc.h */

  /* 
   * Init Inner Hair Cell Model
   */
  /* For Human Only !!!!!!! Hair Cell low pass filter (tdres,cutoff,gain,order) */
  switch(species)
  {
	case 0:
	default: /* Human */
	  initLowPass(&(p->ihc.hclp),tdres,4500.0,1.0,7);
	  break;
	case 1:  /* Cat, this version use the Laurel's revcor data to determine taumin, taumax - only good for low freqs*/
	case 9:  /* Also for universal CAT; as in JASA 2001 Zhang et al. - good for all CFs */
	  initLowPass(&(p->ihc.hclp),tdres,3800.0,1.0,7);
	  break;
  };
  p->ihc.run = runHairCell;
  p->ihc.run2 = runHairCell2;
  p->ihc.hcnl.A0 = 0.1;  /* Inner Hair Cell Nonlinear Function */
  p->ihc.hcnl.B = 2000.;
  p->ihc.hcnl.C = 1.74;
  p->ihc.hcnl.D = 6.87e-9;
  p->ihc.hcnl.run = runIHCNL;
  p->ihc.hcnl.run2 = runIHCNL2;
return;
}

double runAN(TAuditoryNerve *p,double x)
{ 
  double bmout,ihcout,ppiout,sout;
  bmout = p->bm.run(&(p->bm),x);
  ihcout = p->ihc.run(&(p->ihc),bmout);
  ppiout = p->ihcppi.run(&(p->ihcppi),ihcout);
  sout = p->syn.run(&(p->syn),ppiout);
return(sout);  
}

void runAN2(TAuditoryNerve *p,const double *in,double *out,const int length)
{
  p->bm.run2(&(p->bm),in,out,length);
  p->ihc.run2(&(p->ihc),out,out,length);
  p->ihcppi.run2(&(p->ihcppi),out,out,length);
  p->syn.run2(&(p->syn),out,out,length);
return;
}

/* User Get_Tau, initGammaTone, initBoltzman*/
void initBasilarMembrane(TBasilarMembrane* bm,int model, int species, double tdres, double cf)
{ /*
   *
   * ##### Get Basilar Membrane ########
   * 1. Get a structure of BasilarMembrane
   * 2. Specify wide band filter in the BasilarMembrane Model 
   *    //WB filter not used for all model versions, but it's computed here anyway
   * 3. Specify the OHC model in the control path
   *    3.2 Specify the NL function used in the OHC
   * 4. Specify the AfterOHC NL in the control path
   * */
 
  int bmmodel;
  double taumax,taumin,taurange; /* general */
  double x, centerfreq,tauwb,tauwbmin,dtmp,wb_gain; /* for wb */
  double absdb,ohcasym; /* for ohc */
  double dc,R,minR; /* for afterohc */
  TNonLinear *p;

   /* Model numbers as used in ARLO Heinz et al., Fig. 4 */
  if(model == 1) bmmodel = FeedForward_NL;
  else if(model == 2) bmmodel = FeedBack_NL;
  else if(model == 3) bmmodel = Sharp_Linear;
  else if(model == 4) bmmodel = Broad_Linear;
  else if(model == 5) bmmodel = Broad_Linear_High;
  bm->bmmodel = bmmodel;
  bm->tdres = tdres;

  /*
   *  Determine taumax,taumin,order here
   */
  bm->run = runBasilarMembrane;
  bm->run2 = run2BasilarMembrane;

  bm->bmorder = 3;
  Get_tau(species,cf,bm->bmorder,&taumax,&taumin,&taurange);

  bm->TauMax = taumax;
  bm->TauMin = taumin;
  if(bm->bmmodel&Broad_ALL) 
    bm->tau = taumin;
  else
    bm->tau = taumax;
  
  initGammaTone(&(bm->bmfilter),tdres,cf,bm->tau,1.0,bm->bmorder);
  initGammaTone(&(bm->gfagain),tdres,cf,taumin,1.0,1);
  /*
   * Get Wbfilter parameters
   */
  x = cochlea_f2x(species,cf);
  centerfreq = cochlea_x2f(species,x+1.2); /* shift the center freq Qing use 1.1 shift */
  bm->wborder = 3;
  tauwb = taumin+0.2*(taumax-taumin);
  tauwbmin = tauwb/taumax*taumin;
  dtmp = tauwb*TWOPI*(centerfreq-cf);
  wb_gain = pow((1+dtmp*dtmp), bm->wborder/2.0);

  bm->TauWB = tauwb;
  bm->TauWBMin = tauwbmin;
  initGammaTone(&(bm->wbfilter), tdres, centerfreq, tauwb, wb_gain, bm->wborder);
  bm->A=(tauwb/taumax-tauwbmin/taumin)/(taumax-taumin);
  bm->B=(taumin*taumin*tauwb-taumax*taumax*tauwbmin)/(taumax*taumin*(taumax-taumin));

  /*
   * Init OHC model
   */
  absdb = 20; /* The value that the BM starts compression */
  initLowPass(&(bm->ohc.hclp), tdres, 800.0, 1.0, 3); /* This is now use in both Human & Cat MODEL */
  /*/ parameter into boltzman is corner,slope,strength,x0,s0,x1,s1,asym
  // The corner determines the level that BM have compressive nonlinearity */
  ohcasym = 7.0;
  /*/set OHC Nonlinearity as boltzman function combination */
  init_boltzman(&(bm->ohc.hcnl),absdb-12,0.22,0.08,5,12,5,5,ohcasym);
  bm->ohc.run = runHairCell;
  bm->ohc.run2 = runHairCell2;
  
  /*
   * Init AfterOHC model
   */
  p = &(bm->afterohc);
  dc = (ohcasym-1)/(ohcasym+1.0)/2.0;
  dc-=0.05;
  minR = 0.05;
  p->TauMax = taumax;
  p->TauMin = taumin;
  R = taumin/taumax;
  if(R<minR) p->minR = 0.5*R;
  else p->minR   = minR;
  p->A = p->minR/(1-p->minR); /* makes x = 0; output = 1; */
  p->dc = dc;
  R = R-p->minR;
  /* This is for new nonlinearity */
  p->s0 = -dc/log(R/(1-p->minR));
  p->run = runAfterOhcNL;
  p->run2 = runAfterOhcNL2;
return;
}

/** pass the signal through the tuning filter.
    using different model
    Sharp_Linear | Broad_Linear | Broad_Linear_High | FeedBack_NL | FeedForward_NL
*/
double runBasilarMembrane(TBasilarMembrane *bm, double x){
  double out;
  const int length = 1;
  run2BasilarMembrane(bm, &x, &out, length);
return(out);
}

void run2BasilarMembrane(TBasilarMembrane *bm, const double *in, double *out, const int length)
{
  register int i;
  double wb_gain,dtmp,taunow;
  double wbout,ohcout;
  double x, x1,out1;

  for(i=0; i<length; i++)
  {
	  x = in[i];
	  x1 = bm->bmfilter.run(&(bm->bmfilter),x); /*/ pass the signal through the tuning filter */
	  switch(bm->bmmodel){
             default:
             case Sharp_Linear: /* / if linear model, needn't get the control signal */
	     case Broad_Linear:
		out1 = x1; break;
	     case Broad_Linear_High: /* /adjust the gain of the tuning filter */
              	out1 = x1*pow((bm->TauMin/bm->TauMax),bm->bmfilter.Order);
	              break;
	     case FeedBack_NL: /* /get the output of the tuning filter as the control signal */
		wbout = x1;
		
		ohcout = bm->ohc.run(&(bm->ohc),wbout); /*/ pass the control signal through OHC model */
		/*/ pass the control signal through nonliearity after OHC */
		bm->tau = bm->afterohc.run(&(bm->afterohc),ohcout);     
		/*/ set the tau of the tuning filter */
                bm->bmfilter.settau(&(bm->bmfilter),bm->tau);     
		/*  Gain Control of the tuning filter */
              	out1 = pow((bm->tau/bm->TauMax),bm->bmfilter.Order)*x1;  
	              break;
	     case FeedForward_NL:
		/*/get the output of the wide-band pass as the control signal */
		wbout = bm->wbfilter.run(&(bm->wbfilter),x);
		/*/ scale the tau for wide band filter in control path */
		taunow = bm->A*bm->tau*bm->tau-bm->B*bm->tau;

		bm->wbfilter.settau(&(bm->wbfilter),taunow);  /*/set the tau of the wide-band filter*/
		/*/ normalize the gain of the wideband pass filter as 0dB at CF */
		dtmp = taunow*TWOPI*(bm->wbfilter.F_shift-bm->bmfilter.F_shift);
		wb_gain = pow((1+dtmp*dtmp), bm->wbfilter.Order/2.0);
		bm->wbfilter.gain=wb_gain;
		
		ohcout = bm->ohc.run(&(bm->ohc),wbout); /*/ pass the control signal through OHC model*/
	       	/*/ pass the control signal through nonliearity after OHC */
		bm->tau = bm->afterohc.run(&(bm->afterohc),ohcout);
	   	/*/ set the tau of the tuning filter */
	        bm->bmfilter.settau(&(bm->bmfilter),bm->tau);
	      	/*/ Gain Control of the tuning filter */
              	out1 = pow((bm->tau/bm->TauMax),bm->bmfilter.Order)*x1;
              	break;
	  };
	  out[i] = out1;
  };
  bm->gfagain.run2(&(bm->gfagain),out,out,length);

  return;
}


/* ********************************************************************************************
 * * Get TauMax, TauMin for the tuning filter. The TauMax is determined by the bandwidth/Q10
    of the tuning filter at low level. The TauMin is determined by the gain change between high
    and low level
    Also the calculation is different for different species
 */
double Get_tau(int species, double cf,int order, double* _taumax,double* _taumin,double* _taurange)
{
  double taumin,taumax,taurange;
  double ss0,cc0,ss1,cc1;
  double Q10,bw,gain,ratio;
  double xcf, x1000;
  gain = 20+42*log10(cf/1e3);                /*/ estimate compression gain of the filter */
  if(gain>70) gain = 70;
  if(gain<15) gain = 15;
  ratio = pow(10,(-gain/(20.0*order)));      /*/ ratio of TauMin/TauMax according to the gain, order */

  /*/ Calculate the TauMax according to different species */
  switch(species)
    {
    case 1:
      /* Cat parameters: Tau0 vs. CF (from Carney & Yin 88) - only good for low CFs*/
      /* Note: Tau0 is gammatone time constant for 'high' levels (80 dB rms) */
      /* parameters for cat Tau0 vs. CF */
      xcf = cochlea_f2x(species,cf);     /* position of cf unit; from Liberman's map */
      x1000 = cochlea_f2x(species,1000); /* position for 1 Khz */
      ss0 = 6.; cc0 = 1.1; ss1 = 2.2; cc1 = 1.1;
      taumin = ( cc0 * exp( -xcf / ss0) + cc1 * exp( -xcf /ss1) ) * 1e-3;  /* in sec */
      taurange = taumin * xcf/x1000; 
      taumax = taumin+taurange;
      break;
    case 0:
      /* Human Parameters: From Mike Heinz: */
      /* Bandwidths now are based on Glasberg and Moore's (1990) ERB=f(CF,level) equations  */
      taumax =  1./(TWOPI*1.019*erbGM(cf));
      break;
    case 9:
      /* Universal species from data fitting : From Xuedong Zhang,Ian (JASA 2001) */
      /* the Q10 determine the taumax(bandwidths at low level) Based on Cat*/
      Q10 = pow(10,0.4708*log10(cf/1e3)+0.4664);
      bw = cf/Q10;
      taumax = 2.0/(TWOPI*bw);
    }; /*/end of switch */
  taumin =  taumax*ratio;
  taurange = taumax-taumin;
  *_taumin = taumin;
  *_taumax = taumax;
  *_taurange = taurange;
  return 0;
}

/*/ --------------------------------------------------------------------------------
 ** Calculate the location on Basilar Membrane from best frequency
*/
double cochlea_f2x(int species,double f)
{
  double x;
  switch(species)
    {
    case 0: /*/human */
      x=(1.0/0.06)*log10((f/165.4)+0.88);
      break;
    default:
    case 1: /*/cat */
      x = 11.9 * log10(0.80 + f / 456.0);
      break;
    };
   return(x);
}
/* --------------------------------------------------------------------------------
 ** Calculate the best frequency from the location on basilar membrane 
 */
double cochlea_x2f(int species,double x)
{
  double f;
  switch(species)
    {
    case 0: /*/human
      //      if((x>35)||(x<0)) error("BM distance out of human range, [in cochlea_x2f(...)]"); */
      f=165.4*(pow(10,(0.06*x))-0.88);
      break;
    default:
    case 1: /*/cat */
      f = 456.0*(pow(10,x/11.9)-0.80);
      break;
    };
  return(f);
}
/*/ --------------------------------------------------------------------------------
 ** Calculate the erb at 65dB */
double erb51_65(double CF)
{
  double erb1000,erbCf,erb;
  erb1000=24.7*(4.37*1000/1000+1);
  erbCf=24.7*(4.37*CF/1000+1);
  erb=(0.5*erbCf)/(1-(0.38/4000)*erb1000*(65+10*log10(erbCf)-51))+.5*erbCf;
  return(erb);
}
/*/ --------------------------------------------------------------------------------
 ** Calculate the erb using GM's method */
double erbGM(double CF)
{
  double erbCf;
  erbCf=(1/1.2)*24.7*(4.37*CF/1000+1);
  return(erbCf);
}

/*/ ---------------------------------------------------------------------------
 ** Calculate the delay(basilar membrane, synapse for cat*/
double delay_cat(double cf)
{
  /* DELAY THE WAVEFORM (delay buf1, tauf, ihc for display purposes)  */
  /* Note: Latency vs. CF for click responses is available for Cat only (not human) */ 
  /* Use original fit for Tl (latency vs. CF in msec) from Carney & Yin '88 
     and then correct by .75 cycles to go from PEAK delay to ONSET delay */
  double A0 = 8.13; /* from Carney and Yin '88 */
  double A1 = 6.49;
  double x = cochlea_f2x(1,cf); /*/cat mapping */
  double delay = A0 * exp( -x/A1 ) * 1e-3 - 1.0/cf;
return(delay);
}


Loading data, please wait...