Cancelling redundant input in ELL pyramidal cells (Bol et al. 2011)

 Download zip file 
Help downloading and running models
Accession:140465
The paper investigates the property of the electrosensory lateral line lobe (ELL) of the brain of weakly electric fish to cancel predictable stimuli. Electroreceptors on the skin encode all signals in their firing activity, but superficial pyramidal (SP) cells in the ELL that receive this feedforward input do not respond to constant sinusoidal signals. This cancellation putatively occurs using a network of feedback delay lines and burst-induced synaptic plasticity between the delay lines and the SP cell that learns to cancel the redundant input. Biologically, the delay lines are parallel fibres from cerebellar-like granule cells in the eminentia granularis posterior. A model of this network (e.g. electroreceptors, SP cells, delay lines and burst-induced plasticity) was constructed to test whether the current knowledge of how the network operates is sufficient to cancel redundant stimuli.
Reference:
1 . Bol K, Marsat G, Harvey-Girard E, Longtin A, Maler L (2011) Frequency-tuned cerebellar channels and burst-induced LTD lead to the cancellation of redundant sensory inputs. J Neurosci 31:11028-38 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Cerebellum;
Cell Type(s): ELL pyramidal cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program; MATLAB;
Model Concept(s): Dendritic Action Potentials; Bursting; Active Dendrites; Synaptic Plasticity; Long-term Synaptic Plasticity; Learning; Unsupervised Learning; STDP; Biofeedback; Noise Sensitivity;
Implementer(s): Bol, Kieran [kieran_bol at hotmail.com];
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "LIFDAPC.c"
#include "BaysDurhamrand.c"
#define pi 3.141593

int
main(void)
{
    /*---Local variables---*/
    long seed;
    int i,j, GL, numw, nbins, ISIbins;
    double *electroRinput, *lowpassnoisevector, *ISIedges, *PSTHedges, ISIbinwidth,Fs, ISI, *weightsbefore, *weightsafter, numberofT, spmodT;
    
    /* OUTPUTS*/
    double *sptime, *rec, avgw=0, *PSTHhist, *PSTHphase, *ISIhist, firing, burst2,burst4; /* Dynamically allocated output arrays*/
    int spsize, nspikes, n2bursts, n4bursts;
    FILE *outpsth, *outparam, *outISI;
    
    /*---Simulation Parameters */
    double delt=0.01, endtime=1750, freq; /* in seconds */
    int maxj;

    /*---Feedforward parameters---*/
    double tau_m= 0.007, fcut1=500, sigma=0.759, I=0.759*0.759, kappa=0.39;

     /*---Feedback parameters---*/
    double eta=0.0036, eta2, Lambda;
    /*if 2- and 4-spike burst, optimal fit is*/
    double wmax=1.5, tau_w=980, g;
    /*if 2-spike burst only, optimal fit is below, with g=0.87*Lambda; */
    /* double wmax=1, tau_w= 4900,g; */
    
    

    
        
   
    /*---Ask user for frequency---*/
    printf("what is the AM frequency?\n");
    scanf("%lf",&freq);
    getchar();
    
    printf("if global, enter 1\n");
    scanf("%d",&GL);
    getchar();

    /*get seed*/
    printf("enter a seed (negative integer) \n");
    scanf("%d",&seed);
    if(seed>0)
    seed=-1*seed;
    
    getchar();

    numw= floor(1/freq/0.0025);

   

    /*electroreceptor adaptation*/
    if(freq==0.5)
        kappa=0.25;
    if(freq==1)
        kappa=0.27;
    if(freq==2)
        kappa=0.31;

   
   

  


    weightsbefore =(double *) calloc(numw,sizeof(double)); 
    weightsafter =(double *) calloc(numw,sizeof(double)); 
    for(i=0;i<numw;i++)
    {
        weightsbefore[i]=wmax;
    }


    if (GL==1)
    {
        Lambda = 1;
        maxj=3;
    }
    else
    {
        Lambda=0;
        maxj=1;
    }

    g=1.44*Lambda;
    
    /*eta decreasing at high freqs to simulate granule cell burst decay */
/*
    if (freq ==12)
        eta = eta/2;
    if (freq ==16)
        eta=eta/4;
    if (freq==20)
        eta=0.0003;
    if (freq==32)
        eta=0.00015;

    eta2=eta/2;
    */

    /*Lambda decreasing at high freqs model*/
/*   if (freq ==12)  Lambda =0.5; 
    if (freq ==16)  Lambda =0.27;
    if (freq==20)   Lambda =0.05;
    if (freq==32)   Lambda =0;
*/

/*---Initialization*/

    spsize= endtime/delt/tau_m/100;
    sptime=             (double *) calloc(spsize, sizeof(double));
    rec=                (double *) calloc(spsize, sizeof(double));
    electroRinput =     (double *) calloc(endtime/delt/tau_m, sizeof(double));
    lowpassnoisevector= (double *) calloc(endtime/delt/tau_m, sizeof(double));

    Fs=1/delt/tau_m; /* sampling rate in real time*/ 

    
    /*---Simulation loop---*/
    for(j=0;j<maxj;j++)
    { 
        /*This function generates lowpass filtered gaussian noise with cutoff frequency 500Hz, assuming Fs =1/0.01/0.007
        If you change tau_m or delt, you NEED to change the lowpassnoise function*/
        lowpassnoise(endtime/delt/tau_m,seed,lowpassnoisevector);
        
        /*feedforward input*/
        for(i=0;i<endtime/delt/tau_m; i++)
        {
            electroRinput[i]=I+kappa*sin(-2*pi*freq*i*delt*tau_m)+sigma/sqrt(fcut1*2/Fs)*lowpassnoisevector[i];
            if(electroRinput[i]<0)electroRinput[i]=0; /*rectification*/
        }
    
        LIFDAPmodel(electroRinput,endtime/delt/tau_m,freq,weightsbefore,numw,g,Lambda,
        eta,eta2,tau_m,tau_w,wmax,sptime,rec,weightsafter, &nspikes, &n2bursts, &n4bursts); 
   
        for(i=0;i<numw;i++) weightsbefore[i]=weightsafter[i];
        
       

    }
    
    /*Printing average firing and burst rates. If global, also printing average weights*/
    burst4=n4bursts/endtime; /*4-sp burst rate in seconds*/
    burst2=n2bursts/endtime; /*2-sp burst rate in seconds*/
    firing=nspikes/endtime;

    printf("average spike rate is %f\n",firing);
    printf("average 2-spike burst rate is %f\n",burst2);
    printf("average 4-spike burst rate is %f\n",burst4);

    if(GL==1)
    {
        for(i=0;i<numw;i++) avgw += weightsafter[i]/numw;
        printf("The average weight value was %f\n",avgw);
    }
    getchar();

    outparam=fopen("params.txt","w");
        fprintf(outparam,"average spike rate is %f\n",firing);
        fprintf(outparam,"average 2-spike burst rate is %f\n",burst2);
        fprintf(outparam,"average 4-spike burst rate is %f\n",burst4);
        fprintf(outparam,"The average weight value was %f\n",avgw);
     
     fclose(outparam);

 

    
    /*Histogramming firing times. nbins is the number of bins period is divided into*/
    nbins=20; /*number of bins of PSTH over period*/
    ISIbins=50; /*number of ISI bins*/
    ISIbinwidth=4.0; /*Width of each ISI histogram bin in milliseconds*/

    numberofT=endtime*freq; /*number of periods simulated*/
    PSTHedges=(double *) calloc(nbins+1, sizeof(double));
    ISIedges= (double *) calloc(ISIbins+1, sizeof(double));
    PSTHhist=(double *) calloc(nbins, sizeof(double));
    PSTHphase=(double *) calloc(nbins, sizeof(double));
    ISIhist= (double *) calloc(ISIbins, sizeof(double));

    for(i=0; i<nbins;i++)
    {
        PSTHhist[i]=0.0; /*Initializing histogram */
        PSTHedges[i]=i/freq/nbins; /*Dividing period into equal bins for histogramming*/
        PSTHphase[i]=pi/nbins +2*pi/nbins*i;
    }
    PSTHedges[nbins]=1/freq; /*adding last bin*/

    for(i=0;i<ISIbins;i++)
    {
        ISIhist[i]=0;
        ISIedges[i]=ISIbinwidth*i; /*Each bin covers 4 ms */
    }
    ISIedges[ISIbins+1]=ISIbins*ISIbinwidth; /*longest ISI shown in histogram. (200ms with these values)*/

    for(i=0;i<nspikes;i++)
    {
        spmodT=fmod(sptime[i+4]*tau_m,1/freq); /*finding firing time modulo the period (firing times were in units of tau_m)*/
        /*note that sptimes is offset by 4 in LIFDAP so that 4-spike burst searching can start immediately*/

        for(j=0;j<nbins;j++)
        {
            if(PSTHedges[j]<=spmodT && spmodT<PSTHedges[j+1]) PSTHhist[j]+=1.0/numberofT*freq*nbins; /*binning spike times and converting to Hz*/
        } /*when hist has totaled all spikes then multiplying by 1/numberofT will equal average number of spikes per bin*/
            /*To convert to average firing rate in bin, divide by duration of bin = 1/freq/nbins     */

        if(i<nspikes-1)
        {
            ISI=(sptime[i+4+1]-sptime[i+4])*tau_m*1000; /*ISIs are in milliseconds now!*/
            for(j=0;j<ISIbins;j++)
            {
                if(ISIedges[j]<=ISI && ISI<ISIedges[j+1]) ISIhist[j]+=1.0/nspikes/ISIbinwidth; /*Normalized to 1*/ 
            }

        }
    }

    /*print PSTH to file*/
    outpsth=fopen("PSTH.txt","w");
    for(i=0;i<nbins;i++)
    
     fprintf(outpsth, "%f \t %f \n",PSTHphase[i],PSTHhist[i]); /*So PSTHs are always plotted against phase.
                                                                 That way, different frequencies can be easily plotted against each other.*/
    
     fclose(outpsth);

   /*Print ISI to file*/
    outISI=fopen("ISI.txt","w");
    for(i=0;i<ISIbins;i++)
    
     fprintf(outISI, "%f \t %f \n",ISIedges[i]+ISIbinwidth/2,ISIhist[i]);
    
     fclose(outISI);
   
return(0);

}

/*

%histogramming the firing times
npts=20; %20 bins in histogram
sptimes=A(A>0)*tau_m; %all the spike times (A is in units of tau_m)
spper= rem(sptimes,abs(1/freq));  %finding spike times mod period of AM

edges=linspace(0,abs(1/freq),npts+1); 
n=histc(spper,edges);
numper=endtime*freq; %number of periods avgd over to get mean PSTH
n=n/numper*npts*freq; %n/numper = avg #  spikes in the bin in 1 period; 
%T/npts =1/(f*npts) = length of each bin 



%ISI stuff
ISI=sptimes(2:end)-sptimes(1:end-1);
ISI=ISI*1000; %ISI now in ms
ISIn=histc(nearest(ISI*10)/10,0:200); %was 501
ISIn=ISIn/sum(ISIn(1:end-1));


*/