Maximum entropy model to predict spatiotemporal spike patterns (Marre et al. 2009)

 Download zip file 
Help downloading and running models
Accession:125290
This MATLAB code implements a model-based analysis of spike trains. The analysis predicts the occurrence of spatio-temporal patterns of spikes in the data, and is based on a maximum entropy principle by including both spatial and temporal correlations. The approach is applicable to unit recordings from any region of the brain. The code is based on Marre, et al., 2009. The MATLAB code was written by Sami El Boustani and Olivier Marre.
Reference:
1 . Marre O, El Boustani S, Fr├ęgnac Y, Destexhe A (2009) Prediction of spatiotemporal patterns of neural activity from pairwise correlations. Phys Rev Lett 102:138101 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Maximum entropy models;
Implementer(s): El Boustani, Sami [elboustani at unic.cnrs-gif.fr]; Marre, Olivier [marre at unic.cnrs-gif.fr];
%STANDARD PROGRAM FOR DATA ANALYSIS: analyze the file "spikes.txt"

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% General parameters %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

i2mPath;
BinTab=[10 20];%Put here the different bin sizes you want to test
SizeMax=5;
nbpatterns=1000;
n_sample=1000;
N=10;%Enter here the number of neurons to be included in the analysis
PrintFig=1;%Put to 1 to store some figures in the Figures directory

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Initializing storage matrices %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

divjs=zeros(2,length(BinTab),SizeMax,3);
divjsNoJ1=zeros(2,length(BinTab),SizeMax,3);
divjsNoJ=zeros(2,length(BinTab),SizeMax,3);

coefsh= zeros(N,3,length(BinTab));
coefsJ = zeros(N,N,3,length(BinTab));
coefsJ1 = zeros(N,N,3,length(BinTab));
coefshr = zeros(N,3,length(BinTab));
coefsJr = zeros(N,N,3,length(BinTab));
coefshstat = zeros(N,3,length(BinTab));
coefsJstat = zeros(N,N,3,length(BinTab));

PexpTab = zeros(nbpatterns,length(BinTab),SizeMax,3);
PthTab = zeros(nbpatterns,length(BinTab),SizeMax,3);
PthNoJ1Tab = zeros(nbpatterns,length(BinTab),SizeMax,3);
PthNoJTab = zeros(nbpatterns,length(BinTab),SizeMax,3);

r_raw = LoadRaster('spikes.txt');%Enter here the file where the spike times are stored
f=find(r_raw(2,:)<=N);%Reduce the number of neurons to be included in the analysis to N
r = r_raw(:,f);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Loop on the time bin size %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for BinIndex=1:length(BinTab)
    BinSize=BinTab(BinIndex);
    
    % Load the raster data: it is also possible to concatenate the
    % information from several rasters, and to have different types of
    % data (indexed by the variable "state"). 
    raster{1,1} = BinRaster(r,BinSize,0,max(r(1,:)));%First raster for state 1
    raster{2,1} = [];%Second raster for state 1
    raster{3,1} = [];%Third raster for state 1

    % Compute the Correlation matrices and the mean from the rasters, 
    datalenRaw = zeros(3,1);
    for i=1:size(raster,1)
        mAv{i,BinIndex} = zeros(N,1);
        CAv{i,BinIndex} = zeros(N,N);
        C1Av{i,BinIndex} = zeros(N,N);
        for j=1:size(raster,2)
            if ~isempty(raster{i,j})
                datalenRaw(i,j) = size(raster{i,j},2);
                mRaw{i,j} = mean(raster{i,j},2);
                CRaw{i,j} = raster{i,j}*transpose(raster{i,j})/datalenRaw(i,j);
                C1Raw{i,j} = raster{i,j}(:,1:(datalenRaw(i,j)-1))*transpose(raster{i,j}(:,2:datalenRaw(i,j)))/(datalenRaw(i,j)-1);                            
                mAv{i,BinIndex} = mAv{i,BinIndex} + datalenRaw(i,j)*mRaw{i,j};
                CAv{i,BinIndex} = CAv{i,BinIndex} + datalenRaw(i,j)*CRaw{i,j};
                C1Av{i,BinIndex} = C1Av{i,BinIndex} + datalenRaw(i,j)*C1Raw{i,j};
            end
        end
        datalen(i) = sum(datalenRaw(i,:));
        mAv{i,BinIndex} = mAv{i,BinIndex}/datalen(i);
        CAv{i,BinIndex} = CAv{i,BinIndex}/datalen(i);
        C1Av{i,BinIndex} = C1Av{i,BinIndex}/datalen(i);
    end

    % Computing for each state the distribution Pth and Pexp
    for state=1:1%Here there is only one state, but you can have more by just modifying the loading procedure abov. 
        for j=1:size(raster,2)
            rast{j} = raster{state,j};
        end
        m=mAv{state,BinIndex};
        C=CAv{state,BinIndex};
        C1=C1Av{state,BinIndex};
        [mtot,Ctot] = ConvertSpatial(m,C,C1);
        % Tanaka first guess
        [htotguess,Jtotguess]=hJinitial2(mtot,Ctot,2*N);
        [hstatGuess,JstatGuess]=hJinitial2(m,C,N);       
        % Monte-Carlo gradient descent
        [htot,Jtot,time_diff{state,BinIndex}]=GradientMC(mtot,Ctot,htotguess,Jtotguess,1000000,0.1,200,0.005);
        [hstat,Jstat,time_diff_stat{state,BinIndex}]=GradientMC(m,C,hstatGuess,JstatGuess,1000000,0.1,100,0.005);  
        % Storing the corresponding coefficients for the different bin
        % sizes and states
        [h,J,J1,hr,Jr] = DeConvertSpatial(htot,Jtot,hstat,Jstat);
        coefsh(:,state,BinIndex) = h;
        coefsJ(:,:,state,BinIndex) = J;
        coefsJ1(:,:,state,BinIndex) = J1;
        coefshr(:,state,BinIndex) = hr;
        coefsJr(:,:,state,BinIndex) = Jr;
        coefshstat(:,state,BinIndex) = hstat;
        coefsJstat(:,:,state,BinIndex) = Jstat;
    
        % Loop on the template size
        for TempSize=2:SizeMax
            
            %%% Model with spatio-temporal correlations %%%
            [Pexp,Pth,Features,Plist,RasterSize]=OcPred(rast,hstat,Jstat,hr,Jr,h,J,J1,TempSize,nbpatterns);%Pick some patterns, estimate empirically and compare o the prediction of the model. 
            divjs(BinIndex,TempSize,state)=PlotOcFit(Pexp,Pth,Features,BinSize,TempSize,N,RasterSize);% Plot the result
            if (PrintFig==1)
                print('-depsc',[DataFigures,'Data_s',int2str(state), '_b',int2str(BinSize),'_' int2str(N) 'x',int2str(TempSize),'.eps']);
            end
            PexpTab(1:length(Pexp),BinIndex,TempSize,state) = Pexp;
            PthTab(1:length(Pexp),BinIndex,TempSize,state) = Pth;
            PlistTab{BinIndex,TempSize,state} = Plist;
            FeaturesTab(1:length(Pexp),BinIndex,TempSize,state) = Features;
            RasterSizeTab(state,BinIndex) = RasterSize;
            ['Divergence, s' int2str(state) ', b ' int2str(BinSize) ' t ' int2str(TempSize),': ', num2str(divjs(1,BinIndex,TempSize,state)) ]
            close;
            
            %%% Model with spatial correlations %%%
            PthNoJ1 = OcPredFast(Plist,Pexp,hstat,Jstat,zeros(N,1),zeros(N,N),hstat,Jstat,zeros(N,N));
            divjsNoJ1(BinIndex,TempSize,state)=PlotOcFit(Pexp,PthNoJ1,Features,BinSize,TempSize,N,RasterSize);
            if (PrintFig==1)
                print('-depsc',[DataFigures,'DataNoJ1_s',int2str(state), '_b',int2str(BinSize),'_' int2str(N) 'x',int2str(TempSize),'.eps']);
            end
            PthNoJ1Tab(1:length(Pexp),BinIndex,TempSize,state) = PthNoJ1;
            ['Divergence when no Temp, s' int2str(state) ', b ' int2str(BinSize) ' t ' int2str(TempSize),': ', num2str(divjsNoJ1(1,BinIndex,TempSize,state)) ]
            close;
        
            %%% Model without correlation %%%
            PthNoJ = OcPredFast(Plist,Pexp,hstat,zeros(N,N),zeros(N,1),zeros(N,N),hstat,zeros(N,N),zeros(N,N));
            divjsNoJ(BinIndex,TempSize,state)=PlotOcFit(Pexp,PthNoJ,Features,BinSize,TempSize,N,RasterSize);
            if (PrintFig==1)
                print('-depsc',[DataFigures,'DataNoJ_s',int2str(state), '_b',int2str(BinSize),'_' int2str(N) 'x',int2str(TempSize),'.eps']);
            end
            PthNoJTab(1:length(Pexp),BinIndex,TempSize,state) = PthNoJ;
            ['Divergence when no J, s' int2str(state) ', b ' int2str(BinSize) ' t ' int2str(TempSize),': ', num2str(divjsNoJ(1,BinIndex,TempSize,state)) ]
            close;

        end
    end
end

save([pwd,'/Workspace/DataIsing.mat'])

Loading data, please wait...