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];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% BATCH FOR GLAUBER MODEL %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

i2mPath;
SizeMax=5;
nbpatterns=1000;
BinSize = 10;
datalen=100000;
N=8;
state=1;
BinIndex = 1;
tauTab = [1.5 2];

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

PexpTab = zeros(nbpatterns,SizeMax);
PthTab = zeros(nbpatterns,SizeMax);
PthNoJ1Tab = zeros(nbpatterns,SizeMax);
PthNoJTab = zeros(nbpatterns,SizeMax);
PthNoJ0Tab = zeros(nbpatterns,SizeMax);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Generating the Glauber raster and statistics %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%J1test = random('unif',-0.15,0.15,N,N);
J1test = -0.1 + 0.2*rand(N,N);%or: random('unif',-0.1,0.1,N,N);
htest = -1.05 + 0.05*rand(N,1);%or: random('unif',-1.05,-1.0,N,1);

for k=1:length(tauTab)%Here we generate raster with the Glauber model for two different taus
    %Alternatively, a raster can be loaded by using the LoadRaster and
    %BinRaster functions
    tau = tauTab(k);
    raster{k} = GlauberRaster(J1test,htest,N,datalen,1,tau);%For this parameter tau, we generate a raster
    
    %HERE START THE ANALYSIS WITH THE STATISTICAL MODEL. In this example, the data are in
    %raster{k}, which is a binned raster composed of -1 and 1. See the
    %sub-directory /LoadRaster to load from data instead of using the
    %Glauber model. 
    m{k} = mean(raster{k},2);
    C{k} = raster{k}*transpose(raster{k})/datalen;
    C1{k} = raster{k}(:,1:(datalen-1))*transpose(raster{k}(:,2:datalen))/(datalen-1);
    [mtot,Ctot] = ConvertSpatial(m{k},C{k},C1{k});
    % First guess with the Tanaka approximation
    [htotguess,Jtotguess]=hJinitial2(mtot,Ctot,2*N); 
    [hstatGuess,JstatGuess]=hJinitial2(m{k},C{k},N);       
    
    % Monte Carlo gradient descent
    [htot,Jtot]=GradientMC(mtot,Ctot,htotguess,Jtotguess,10000000,0.1,50,0.007); 
    [hstat,Jstat]=GradientMC(m{k},C{k},hstatGuess,JstatGuess,1000000,0.1,100,0.007);
    [h,J,J1,hr,Jr] = DeConvertSpatial(htot,Jtot,hstat,Jstat); 
    
    %Store the fitted coefficients h and J in some tables for further use. 
    coefsh(:,k,BinIndex) = h;
    coefsJ(:,:,k,BinIndex) = J;
    coefsJ1(:,:,k,BinIndex) = J1;
    coefshr(:,k,BinIndex) = hr;
    coefsJr(:,:,k,BinIndex) = Jr;
    coefshstat(:,k,BinIndex) = hstat;
	coefsJstat(:,:,k,BinIndex) = Jstat;
        
    %Estimating the 3 models performance: predicting probabilities of
    %individual patterns (figure 1)
    for TempSize=1:SizeMax%For patterns of size 1 to SizeMax
    
        %%% Model with spatio-temporal correlations %%%
        rast{1} = raster{k};
        [Pexp,Pth,Features,Plist,RasterSize]=OcPred(rast,hstat,Jstat,hr,Jr,h,J,J1,TempSize,nbpatterns);%See OcPred picks up some patterns, compute their ocurrence empirically, and predict it from the model. 
    	divjs(TempSize,k)=PlotOcFit(Pexp,Pth,Features,BinSize,TempSize,N);% Plot the prediction vs empirical estimations
        
        PexpTab(1:length(Pexp),BinIndex,TempSize,k) = Pexp;% These 3 lines store the results in table for subsequent use in Fig 1.
        PthTab(1:length(Pexp),BinIndex,TempSize,k) = Pth;
        PlistTab{BinIndex,TempSize,k} = Plist;
        FeaturesTab(1:length(Pexp),BinIndex,TempSize,k) = Features;
        ['Divergence, t ' int2str(TempSize),': ', num2str(divjs(1,TempSize,k)) ]
    	close;

        %%% Model with only spatial correlations %%%    
        PthNoJ1 = OcPredFast(Plist,Pexp,hstat,Jstat,zeros(N,1),zeros(N,N),hstat,Jstat,zeros(N,N));
    	divjsNoJ1(TempSize,k)=PlotOcFit(Pexp,PthNoJ1,Features,BinSize,TempSize,N);
    	PthNoJ1Tab(1:length(Pexp),BinIndex,TempSize,k) = PthNoJ1;
        ['Divergence when no Temp, t ' int2str(TempSize),': ', num2str(divjsNoJ1(1,TempSize,k)) ]
        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(TempSize,k)=PlotOcFit(Pexp,PthNoJ,Features,BinSize,TempSize,N);
        PthNoJTab(1:length(Pexp),BinIndex,TempSize,k) = PthNoJ;
        ['Divergence when no J, t ' int2str(TempSize),': ', num2str(divjsNoJ(1,TempSize,k)) ]
        close;
    
    end

    
end


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

Loading data, please wait...