Boris Gourevich, Jos J. Eggermont

Maximum decoding abilities of temporal patterns and synchronized firings: application to auditory neurons responding to click trains and amplitude modulated white noise

Boris Gourevich, Jos J. Eggermont

Department of Physiology and Biophysics, Department of Psychology,

University of Calgary, Calgary, Alberta, Canada.

Supplementary material

A Matlab 7 code is provided for the extraction of patterns (extract_patterns) and the selection of variables based on the Forward Method and the Mutual Information (MI_maximization). Two short programs show how to use extract_patterns and MI_maximization. Alternative versions for patterns extraction (based on binning the time axis) and MI_maximization (based on Backward selection of variables) are also provided at the end.

Main programs

Extract_patterns

function [patterns_ch_archive,patterns_seq_archive,patterns_histo_archive,patterns_ind_sub_archive,...

patterns_special_value]=extract_patterns(event_times,binw,d1,d2,d3,option)

%extract all temporal patterns from spike events of a set of channels

%version binning the time difference between events

%

%USE :

%[patterns_ch_archive,patterns_seq_archive,patterns_histo_archive,patterns_

% ind_sub_archive,patterns_special_value]=extract_patterns(event_times,binw,d1,d2,d3)

%

%INPUT : event_times - matrice p*3 with spike time (1st column),

% channels (2nd column) and spike times trigerred by the stimulus (3rd column)

% binw - bin width

% d1 - minimum number of repeats for a pattern to be kept (d1 of rule 1 - see paper)

% d2 - maximum duration (in ms) of a pattern to be kept (d2 of rule 2 - see paper)

% d3 - maximum number of appearances of a chanel in a pattern (d3 of rule 3 - see paper)

% option - (optional) values needed to analyze the PSTH (Modulation Frequency in the paper)

%

%OUTPUT : patterns_ch_archive - channels involved in patterns. Cell array (one cell per pattern size).

% An array associated with patterns of size n contains the list of n-uplets of channels

% involved in each pattern

% patterns_seq_archive - delays between channels within patterns (in s). Cell array (one cell

% per pattern size).

% List of delays of channels (relative to the previous one) corresponding with the patterns

% in patterns_ch_archive

% patterns_histo_archive - occurences of patterns corresponding with those in patterns_ch_archive.

% Cell array (one cell per pattern size).

% patterns_ind_sub_archive - indices of patterns in patterns_ch_archive{s-1} having generated the

% patterns in patterns_ch_archive{s-1}

% Cell array (one cell per pattern size).

% patterns_special_value - value of each pattern for a given

% parameter (if input "option" is used). Synchronization Rate was

% used in the paper. Cell array (one cell per pattern size).

%% Checking input

if nargin<6,

option=0;

end;

if nargin<5,

error('Not enough input arguments. See help extract_patterns');

end;

%% General parameters

tic

%maximum number of channels

n_ch=max(event_times(:,2));

%convert d2 into number of bins

d2=floor(d2/binw);

%size of buffer - should be ok on most PCs but have to be lowered otherwise

taille=200000;%number of patterns

taille2=5000000;%number of patterns * occurences of each one

%FR of each channel

fr=hist(event_times(:,2),1:n_ch);

%select channels

ch_involved=find(fr>d1);

%% Patterns of size 2

% initialization of matrices for patterns of size 2

index_new=1:size(event_times,1)-1;

base_pat=zeros(taille2,3,'uint8');

%matrices of information related to each occurence of patterns

patterns_time_current_special=zeros(taille2,1,'single');

patterns_index_current=zeros(taille2,1,'uint32');

patterns_index_current_end=zeros(taille2,1,'uint16');

l_event_times=size(event_times,1);

%loop: each index N of event_times is paired with the index N+compt. All

%index are taken into account simultaneously through the set of index

%index_new

compt=0;

no_stop=1;

count_tot_patterns=1;

while no_stop,

compt=compt+1;

%prevent using index exceeding matrix dimension

if index_new(end)+compt>l_event_times,

index_new(end)=[];

end;

%time delay between the two events

base_pat_temp=event_times(index_new+compt,1)-event_times(index_new,1);

%delay is binned

base_pat_temp=uint8(floor(base_pat_temp/binw));

%apply rule 2

ind=logical(base_pat_temp<=d2);

%if the index N+compt generates a delay >d2 with index N, then so will

%do the index n+compt+1. Such cases can be removed for the next loop.

%Remaining index are those to store for this current loop.

index_new=index_new(ind);

nb=length(index_new);

if nb>0,

%store pairs

base_pat(count_tot_patterns:count_tot_patterns+nb-1,:)=[base_pat_temp(ind) ...

event_times(index_new,2) event_times(index_new+compt,2)];

%index of first event of pair for each pair

patterns_index_current(count_tot_patterns:count_tot_patterns+nb-1,:)=index_new;

%index difference between the last event and the first event of

%the pattern

patterns_index_current_end(count_tot_patterns:count_tot_patterns+nb-1,:)=compt;

%stimulus trigered time

patterns_time_current_special(count_tot_patterns:count_tot_patterns+nb-1,:)=event_times(index_new,3);

count_tot_patterns=count_tot_patterns+nb;

else

no_stop=0;

end;

end;

%remove the useless zeros

base_pat=base_pat(1:count_tot_patterns-1,:);

patterns_time_current_special=patterns_time_current_special(1:count_tot_patterns-1)';

patterns_index_current=patterns_index_current(1:count_tot_patterns-1)';

patterns_index_current_end=patterns_index_current_end(1:count_tot_patterns-1)';

%find repeated patterns among pairs

[patterns_ch_current,b,index_unique_patterns]=unique(base_pat,'rows');%store channel information

patterns_ch_current=patterns_ch_current';

nb_current=size(patterns_ch_current,2);

patterns_seq_current=[zeros(1,nb_current,1);patterns_ch_current(1,:)];%store delay information

patterns_ch_current(1,:)=[];

patterns_histo_current=hist(index_unique_patterns,1:nb_current);%store count information

%apply rule 1

ind_to_keep=find(patterns_histo_current>d1);

patterns_ch_current=patterns_ch_current(:,ind_to_keep);

patterns_seq_current=patterns_seq_current(:,ind_to_keep);

patterns_histo_current=patterns_histo_current(ind_to_keep);

%keep information about each occurence for selected patterns only

%optimized version using advanced matrice manipulations under Matlab

[ind_to_keep2,index_unique_patterns]=ismember(index_unique_patterns,ind_to_keep);

patterns_index_current=patterns_index_current(ind_to_keep2);

patterns_index_current_end=patterns_index_current_end(ind_to_keep2);

patterns_time_current_special=patterns_time_current_special(ind_to_keep2);

index_unique_patterns=index_unique_patterns(ind_to_keep2);

[temp,index]=sortrows([index_unique_patterns patterns_index_current']);

index_unique_patterns=temp(:,1)';

patterns_index_current=temp(:,2)';

patterns_index_current_end=patterns_index_current_end(index);

patterns_time_current_special=patterns_time_current_special(index);

clear temp;

%first index in patterns_index_current for the n-th pattern

patterns_ind_occ_current=[1 1+cumsum(patterns_histo_current(1:end-1))];

%The basic principle is to get the patterns of size s+1 by aggregating a pattern of size s to a pattern

%of size 2: we now store the information of size 2 to use it as a reference

patterns_ch_ref=patterns_ch_current;

patterns_seq_ref=patterns_seq_current;

patterns_histo_ref=patterns_histo_current;

patterns_index_ref=patterns_index_current;

patterns_index_ref_end=patterns_index_current_end;

patterns_ind_occ_ref=patterns_ind_occ_current;

clear patterns_ind_occ_current;

patterns_ind_pat_ref=index_unique_patterns;

%number of patterns of reference

nb_ref=length(patterns_histo_ref);

%fill in the output variables

clear patterns_seq_archive

%patterns of size 1 (single channels firing once)

fprintf('(Size) : (Number of patterns)\n')

fprintf(['1 : ' num2str(length(ch_involved)) '\n'])

patterns_ch_archive{1}=ch_involved;

patterns_histo_archive{1}=fr(fr>d1);

%patterns of size 2

patterns_ch_archive{2}=patterns_ch_current;

patterns_seq_archive{2}=single(diff(patterns_seq_current))*binw;

patterns_histo_archive{2}=patterns_histo_current;

[a,patterns_ind_sub_archive{2}]=ismember(patterns_ch_current(1,:),patterns_ch_archive{1});

fprintf(['2 : ' num2str(length(patterns_histo_current)) '\n'])

if option>0

%We compute the Synchronization rate from the PSTH (Discrete Fourier

%Transform of the PSTH at the Modulation Frequency), i.e. the stimulus

%trigerred time information stored in patterns_time_current_special

%Compute here your parameter extracted from PSTH

patterns_special_value{1}=abs(accumarray(event_times(:,2),...

exp(-2*i*pi*event_times(:,3)'*option)));

patterns_special_value{1}=patterns_special_value{1}(fr>d1);

patterns_special_value{2}=abs(accumarray(patterns_ind_pat_ref',...

exp(-2*i*pi*patterns_time_current_special*option)));

clear patterns_time_current_special;

else

patterns_special_value=0;

end;

%% Patterns of size > 2

non_stop=1;

%let's start at size 3

s=2;

while non_stop,%still finding some patterns

s=s+1;

%initializations for the current size

patterns_index_new=zeros(1,taille2,'uint32');

patterns_index_end_new=zeros(1,taille2,'uint16');

patterns_index_ind_new=zeros(1,taille2,'uint32');

patterns_ch_seq_new=zeros(s,taille,'uint8');

patterns_seq_new=zeros(s,taille,'uint8');

patterns_histo_new=zeros(1,taille,'uint32');

patterns_ind_sub_new=zeros(1,taille,'uint32');

taille_index=taille2;

count_patterns=1;

count_tot_patterns=1;

%loop on the patterns of reference (size 2)

for j=1:nb_ref,

%variables for the current pattern of reference

ch_ref_current=patterns_ch_ref(:,j);

seq_ref_current=patterns_seq_ref(:,j);

ind_occ_ref_current=patterns_ind_occ_ref(j);

%search potential patterns of size s+1 among patterns of size s: the first channel of the pattern

%of reference must be the last one of the pattern of size s AND the duration of the potential new

%pattern must not be > d2 AND the pattern of size s+1 must not already contain d3 times the same

%channel (equivalent to the pattern of size s must not already contain d3-1 times the last channel

%of the pattern of reference)

ind_sub_new=find(ch_ref_current(1)==patterns_ch_current(end,:)&patterns_seq_current(end,:)+...

seq_ref_current(end)<=d2&sum(patterns_ch_current==ch_ref_current(2))<d3);

%occurences index of the current pattern of reference

index_ref=patterns_index_ref(ind_occ_ref_current:ind_occ_ref_current+...

patterns_histo_ref(j)-1);

%index difference between the first and the last events of the pattern index_ref_end=patterns_index_ref_end(ind_occ_ref_current:ind_occ_ref_current+...

patterns_histo_ref(j)-1);

%Optimization: the following instructions allows managing all patterns of size

%s+1 simultaneously

%get the occurences index of patterns of size s previously selected

temp_ind=find(ismember(index_unique_patterns,ind_sub_new));

patterns_index_current_temp=patterns_index_current(temp_ind);

patterns_index_current_end_temp=patterns_index_current_end(temp_ind);

%which pattern do these occurence index refer to

index_unique_patterns_temp=index_unique_patterns(temp_ind);

%select occurences which match those of the pattern of

%reference, i.e. which will be the occurences of the patterns of

%size s+1: intersection between the index of the last event of the

%patterns and the index of the first event of the pattern of

%reference

[ind_intersect,value_intersect]=ismember(patterns_index_current_temp+...

uint32(patterns_index_current_end_temp),index_ref);

value_intersect=uint32(value_intersect);

%compute the number of occurences for all potential patterns of

%size s+1

histo_new=accumarray(index_unique_patterns_temp',ind_intersect);

histo_new=histo_new(ind_sub_new);

%Apply rule 1

ind=histo_new>d1;

histo_new=histo_new(ind);

ind_sub_new=ind_sub_new(ind);

if ~isempty(ind_sub_new),

%number of patterns of size s+1

nb=length(ind_sub_new);

%store the information for the patterns of size s+1

patterns_ch_seq_new(:,count_patterns:count_patterns+nb-1)=...

[patterns_ch_current(:,ind_sub_new);ch_ref_current(2)*ones(1,nb,'uint8')];

patterns_seq_new(:,count_patterns:count_patterns+nb-1)=[patterns_seq_current(:,ind_sub_new);...

seq_ref_current(end)+patterns_seq_current(end,ind_sub_new)];

patterns_histo_new(count_patterns:count_patterns+nb-1)=histo_new';

patterns_ind_sub_new(count_patterns:count_patterns+nb-1)=ind_sub_new;

%check if storing the new patterns will overflow the current

%matrices. If so, enlarge them.

if count_tot_patterns+sum(histo_new)-1>taille_index,

patterns_index_new=[patterns_index_new zeros(1,taille2,'uint32')];

patterns_index_ind_new=[patterns_index_ind_new zeros(1,taille2,'uint32')];

patterns_index_end_new=[patterns_index_end_new zeros(1,taille2,'uint16')];

taille_index=taille_index+taille2;

end;

ind_intersect=ind_intersect&ismember(index_unique_patterns_temp,ind_sub_new);

patterns_index_new(count_tot_patterns:count_tot_patterns+sum(histo_new)-1)=...

patterns_index_current_temp(ind_intersect);

patterns_index_ind_new(count_tot_patterns:count_tot_patterns+sum(histo_new)-1)=...

count_patterns-1+cumsum([1 diff(index_unique_patterns_temp(ind_intersect))>0]);

patterns_index_end_new(count_tot_patterns:count_tot_patterns+sum(histo_new)-1)=...

patterns_index_current_end_temp(ind_intersect)+index_ref_end(value_intersect(ind_intersect));

count_tot_patterns=count_tot_patterns+sum(histo_new);

count_patterns=count_patterns+nb;

end;

end;

%if some patterns of size s+1 were found

if count_patterns-1>0,

%archive results

clear patterns_index_currentpatterns_index_current_endindex_unique_patterns

patterns_ch_archive{s}=patterns_ch_seq_new(:,1:count_patterns-1);

patterns_seq_archive{s}=patterns_seq_new(:,1:count_patterns-1);

patterns_histo_archive{s}=patterns_histo_new(1:count_patterns-1);

patterns_ind_sub_archive{s}=patterns_ind_sub_new(1:count_patterns-1);

fprintf([num2str(s) ' : ' num2str(count_patterns-1) '\n'])

clear patterns_ch_seq_newpatterns_seq_newpatterns_histo_newpatterns_ind_sub_new

%the information for the patterns of size s allowing to find patterns of size s+1 is replaced by

%that for the patterns of size s+1

patterns_ch_current=patterns_ch_archive{s};

patterns_seq_current=patterns_seq_archive{s};

patterns_seq_archive{s}=single(diff(patterns_seq_archive{s}))*binw;

patterns_index_current=patterns_index_new(1:count_tot_patterns-1);

clear patterns_index_new

patterns_index_current_end=patterns_index_end_new(1:count_tot_patterns-1);

clear patterns_index_end_new

index_unique_patterns=patterns_index_ind_new(1:count_tot_patterns-1);

clear patterns_index_ind_new

if option>0

patterns_special_value{s}=abs(accumarray(index_unique_patterns',...

exp(-2*i*pi*event_times(patterns_index_current,3)*option)));

end;

else

%no patterns found, stop the loop

non_stop=0;

end;

end;

toc

Test of extract_patterns

%% Simulate PSTHs for n channels in response to clicks at 10Hz

ch_number=16;

%number of repetitions

repeats=50;

%click train frequency

frequency=10;

%firing rate for each channel

FR=[50 50 50 50 25 25 25 10 10 10 10 10 10 10 10 2];

data=[];

for i=1:ch_number,

%spikes every 0.1s selected randomly

temp=1/frequency*floor(repeats*FR(i)*rand(repeats*FR(i),1));

%times following an exponential distribution are added to each spike

temp=temp+exprnd(1/(3*FR(i)),repeats*FR(i),1);

%(channels with high FR are made bursting more)

%channel information

data=[data;temp i*ones(repeats*FR(i),1)];

end;

%compute the PSTH

data(:,3)=mod(data(:,1),1);

%sort times

data=sortrows(data,1);

%display PSTH for all channels taken into account

figure;

hist(data(:,3),[0:0.01:1])

%% Test extract_patterns

%rules

d1=10;

d2=0.01;

d3=4;

%bin width (in s)

binw=0.002;

[patterns_ch_archive,patterns_seq_archive,patterns_histo_archive,patterns_ind_sub_archive,...

patterns_special_value]=extract_patterns(data,binw,d1,d2,d3,frequency);

MI maximization by Forward Method

function [MI_est,select_col,history_MI_est]=MI_maximization(distrib,base)

%maximum of mutual information available in a matrix by selecting the

%columns from a forward method

%

%USE :

%[MI_est,select_col,history_MI_est]=MI_maximization(distrib,base)

%

%INPUT : distrib - joint distribution between stimuli in rows and

% observations in columns

% base - log base

%

%OUTPUT : MI_est - final MI estimate

% select_col - selected columns

% history_MI_est - history of optimization: MI and index of

% variable added at each step

%

%Note: distrib can have null values

%% Checking input

if isempty(distrib),

error('Empty distribution');

end;

%% Initialize and normalize some parameters

[ncellx,ncelly]=size(distrib);

history_MI_est=zeros(2,ncelly);

%algorithm stops if the increase of MI is smaller than tolerance: can be

%put to 0 for real data. For simulated MTFs, it generally avoids to add

%some redundant MTFs that would keep MI stable

tolerance=5*10^-5;

%give each column the same probability -> margin for dimension y is one

distrib=distrib./repmat(sum(distrib,1),ncellx,1);

%% Initialization of the Forward selection: choose variable of min entropy

%estimate entropy for each column

H_i=-distrib.*log(distrib);

%deal with NaN values stemming from null values in distrib

H_i=nansum(H_i);

%look for the minimum

[H_i_min,ind_H_i_min]=min(H_i);

%initialize the history of MI maximization

history_MI_est(:,1)=[0;ind_H_i_min];

MI_est=0;

%% Prepare the loop on other columns

%remaining columns

ind_rest=setdiff(1:ncelly,ind_H_i_min);

%Optimization: we decompose I(X,Y)=H(X)+H(Y)-H(X,Y) and

%compute everything related to margin of y first. Compared to the

%pseudo-code, we treat all columns simultaneously (no loop FOR..END)

%-H(X,Y): sum only on x of p(x,y).log(p(x,y)) is equivalent to -H(X,Y)

%before integration over y

mat_ref_y=H_i;

clear H_i;

%Since margin for y is one, this is actually equivalent to

%H(Y)-H(X,Y) before integration over y

%margin_x of the matrix of selected MTFs

margin_x=distrib(:,ind_H_i_min);

%% Loop while MI continues to increase

stop=0;

counter=1;

while ~stop,

%number of steps = normalization for the distribution composed of

%selected variables plus one given that the margins for y are one in

%distrib

counter=counter+1;

%compute H(Y)-H(X,Y) for each column added (sum on previously selected

%columns plus what was in each y value)

mat_ref_y_sum=sum(mat_ref_y(history_MI_est(2,1:counter-1)))+...

mat_ref_y(ind_rest);

mat_ref_y_sum=-mat_ref_y_sum./counter;

%add one column to margin_x

margin_x_new=repmat(margin_x,1,length(ind_rest))+distrib(:,ind_rest);

%normalization of the distributions

margin_x_new=margin_x_new./repmat(sum(margin_x_new),ncellx,1);

%MI after adding i-th column

%compute the new H(X) then add H(Y)-H(X,Y)

MI_i=-margin_x_new.*log(margin_x_new);

MI_i=nansum(MI_i)+mat_ref_y_sum;

%MI in some base

MI_i=MI_i/log(base);

%look for the maximum

[MI_est_max,ind_MI_est_max]=max(MI_i);

if MI_est_max<MI_est+tolerance,

%the adding of any variable leads to a decrease of the MI: the

%algorithm stops

stop=1;

history_MI_est(:,counter)=[MI_est_max;0];

else

%the reference value of MI is updated

MI_est=MI_est_max;

%the variable removed is stored in the history

history_MI_est(:,counter)=[MI_est;ind_rest(ind_MI_est_max)];

%margin_x is updated by adding the distribution of the chosen

%variable

margin_x=margin_x+distrib(:,ind_rest(ind_MI_est_max));

%vector of remaining variables updated

ind_rest(ind_MI_est_max)=[];

%if no variable left, the algorithm stops

if isempty(ind_rest),

stop=1;

end;

end;

end;

%selected columns

if isempty(ind_rest),

select_col=sort(history_MI_est(2,1:counter));

else

select_col=sort(history_MI_est(2,1:counter-1));

end;

history_MI_est=history_MI_est(:,1:counter);

Test of MI_maximization

%% Simulate distributions of PSTHs

%number of observations (patterns in the paper) and stimuli (modulation

%frequencies)

nb_obs=200;

nb_stim=20;

%generation of nb_obs white gaussian noises

data=randn(nb_obs,nb_stim);

%lowpass filtering

data=filter(ones(1,5),1,data,[],2);

data=data(:,end:-1:1);

%normalize each TMTF

data=(data-repmat(min(data,[],2),1,nb_stim)+1);%keep them >0

data=data./repmat(sum(data,2),1,nb_stim);

%three TMTFs

figure;plot(data(25:27,:)')

%% selection of variables

[info_max,ind_info_max,evol_info_max]=MI_maximization(data',nb_stim);

%display optimal TMTFs

figure;plot(data(ind_info_max,:)')

Alternative versions

Extract_patterns2

function [patterns_ch_archive,patterns_seq_archive,patterns_histo_archive,patterns_ind_sub_archive,...

patterns_special_value]=extract_patterns2(event_times,binw,d1,d2,d3,option)

%extract all temporal patterns from spike events of a set of channels

%version binning the time axis (see discussion in the paper)

%

%USE :

%[patterns_ch_archive,patterns_seq_archive,patterns_histo_archive,patterns_

% ind_sub_archive,patterns_special_value]=extract_patterns2(event_times,binw,d1,d2,d3)

%

%INPUT : event_times - matrice p*3 with spike time (1st column),

% channels (2nd column) and spike times trigerred by the stimulus (3rd column)

% binw - bin width

% d1 - minimum number of repeats for a pattern to be kept (d1 of rule 1 - see paper)

% d2 - maximum duration (in ms) of a pattern to be kept (d2 of rule 2 - see paper)

% d3 - maximum number of appearances of a chanel in a pattern (d3 of rule 3 - see paper)

% option - (optional) values needed to analyze the PSTH (Modulation Frequency in the paper)

%

%OUTPUT : patterns_ch_archive - channels involved in patterns. Cell array (one cell per pattern size).

% An array associated with patterns of size n contains the list of n-uplets of channels

% involved in each pattern

% patterns_seq_archive - delays between channels within patterns (in bins). Cell array (one cell

% per pattern size).

% List of delays of channels (relative to the first one) corresponding with the patterns

% in patterns_ch_archive

% patterns_histo_archive - occurences of patterns corresponding with those in patterns_ch_archive.

% Cell array (one cell per pattern size).

% patterns_ind_sub_archive - indices of patterns in patterns_ch_archive{s-1} having generated the

% patterns in patterns_ch_archive{s-1}

% Cell array (one cell per pattern size).

% patterns_special_value - value of each pattern for a given

% parameter (if input "option" is used). Synchronization Rate was

% used in the paper. Cell array (one cell per pattern size).

%% Checking input

if nargin<6,

option=0;

end;

if nargin<5,

error('Not enough input arguments. See help extract_patterns');

end;

%% General parameters

tic

%maximum number of channels

n_ch=max(event_times(:,2));

%convert d2 into number of bins

d2=floor(d2/binw);

%size of buffer - should be ok on most PCs but have to be lowered otherwise

taille=200000;%number of patterns

taille2=5000000;%number of patterns * occurences of each one

%binning of spike times

event_times_binned=[floor(event_times(:,1)/binw) event_times(:,2) floor(event_times(:,3)/binw)];

%fix a rounding impreciseness of matlab

temp_var=event_times(:,1)/binw-round(event_times(:,1)/binw);

ind=find(temp_var<0&temp_var>-10^-4);

temp_var=event_times(:,3)/binw-round(event_times(:,3)/binw);

ind2=find(temp_var<0&temp_var>-10^-4);

clear temp_var

event_times_binned(ind,1)=event_times_binned(ind,1)+1;

event_times_binned(ind2,3)=event_times_binned(ind2,3)+1;

event_times_binned=unique(event_times_binned,'rows');

only_times_binned=unique(event_times_binned(:,[1 3]),'rows');

%FR of each channel

fr=hist(event_times_binned(:,2),1:n_ch);

%select channels

ch_involved=find(fr>d1);

%% Patterns of size 2

% initialization of matrices for patterns of size 2

patterns_ch_current=zeros(2,taille,'uint8');

patterns_seq_current=zeros(2,taille,'uint8');

patterns_histo_current=zeros(1,taille,'uint32');

patterns_ind_occ_current=zeros(1,taille,'uint32');

patterns_time_current=zeros(1,taille2,'uint32');

patterns_ind_pat_current=zeros(1,taille2,'uint32');

patterns_time_current_special=zeros(1,taille2,'single');

%extract patterns of size 2

count_patterns=0;

count_tot_patterns=1;

for ch1=ch_involved,

%select times for channel ch1

times_ch1=event_times_binned(event_times_binned(:,2)==ch1,[1 3]);

for ch2=ch_involved,

%select times for channel ch2

times_ch2=event_times_binned(event_times_binned(:,2)==ch2,[1 3]);

%for each pattern size

for p_size=0:d2,

if p_size~=0||ch1<ch2%do not count things twice

%delay times_ch1 and find corresponding times for ch2

[a,b]=intersect(times_ch1(:,1)+p_size,times_ch2(:,1));

nb=length(a);

%if enough patterns

if nb>d1,

count_patterns=count_patterns+1;%index of this pattern

patterns_ch_current(:,count_patterns)=[ch1;ch2];%channels involved

patterns_seq_current(:,count_patterns)=[0;p_size];%delay

patterns_histo_current(count_patterns)=nb;%number of occurences

%index in patterns_time_current where all occurences of

%this pattern are listed

patterns_ind_occ_current(count_patterns)=count_tot_patterns;

%time of each occurence of this pattern

patterns_time_current(count_tot_patterns:count_tot_patterns+nb-1)=a;

%index of this pattern to which each occurence refers to

patterns_ind_pat_current(count_tot_patterns:count_tot_patterns+nb-1)=count_patterns;

if option>0

patterns_time_current_special(count_tot_patterns:count_tot_patterns+nb-1)...

=times_ch1(b,2)+p_size;%time related to the beginning of the stimulus (PSTH)

end;

count_tot_patterns=count_tot_patterns+nb;

end;

end;

end;

end;

end;

%remove the useless zeros

patterns_ch_current=patterns_ch_current(:,1:count_patterns);

patterns_seq_current=patterns_seq_current(:,1:count_patterns);

patterns_histo_current=patterns_histo_current(1:count_patterns);

patterns_time_current=patterns_time_current(1:count_tot_patterns-1);

if option>0

patterns_time_current_special=patterns_time_current_special(1:count_tot_patterns-1);

end;

patterns_ind_pat_current=patterns_ind_pat_current(1:count_tot_patterns-1);

%The basic principle is to get the patetrns of size s+1 by aggregating a pattern of size s to a pattern

%of size 2: we now store the information of size 2 to use it as a reference

patterns_ch_ref=patterns_ch_current;

patterns_seq_ref=patterns_seq_current;

patterns_histo_ref=patterns_histo_current;

patterns_time_ref=patterns_time_current;

patterns_ind_occ_ref=patterns_ind_occ_current(1:count_patterns);

clear patterns_ind_occ_current;

patterns_ind_pat_ref=patterns_ind_pat_current;

%number of patterns of reference

nb_ref=length(patterns_histo_ref);

%fill in the output variables

clear patterns_seq_archive

%patterns of size 1 (single channels firing once)

fprintf('(Size) : (Number of patterns)\n')

fprintf(['1 : ' num2str(length(ch_involved)) '\n'])

patterns_ch_archive{1}=ch_involved;

patterns_histo_archive{1}=fr(fr>d1);

%patterns of size 2

patterns_ch_archive{2}=patterns_ch_current;

patterns_seq_archive{2}=patterns_seq_current;

patterns_histo_archive{2}=patterns_histo_current;

patterns_ind_sub_archive{2}=patterns_ch_current(1,:);

fprintf(['2 : ' num2str(length(patterns_histo_current)) '\n'])

if option>0

%We compute the Synchronization rate from the PSTH (Discrete Fourier

%Transform of the PSTH at the Modulation Frequency)

%Compute here your parameter extracted from PSTH

patterns_special_value{2}=abs(accumarray(patterns_ind_pat_ref',...

exp(-2*i*pi*patterns_time_current_special*option*binw)));

clear patterns_time_current_special;

else

patterns_special_value=0;

end;

%% Patterns of size > 2

non_stop=1;

%let's start at size 3

s=2;

while non_stop,%still finding some patterns

s=s+1;

%initializations for the current size

patterns_time_new=zeros(1,taille2,'uint32');

patterns_time_ind_new=zeros(1,taille2,'uint32');

patterns_ch_seq_new=zeros(s,taille,'uint8');

patterns_seq_new=zeros(s,taille,'uint8');

patterns_histo_new=zeros(1,taille,'uint32');

patterns_ind_sub_new=zeros(1,taille,'uint32');

taille_time=taille2;

count_patterns=1;

count_tot_patterns=1;

%loop on the patterns of reference (size 2)

for j=1:nb_ref,

%variables for the current pattern of reference

ch_ref_current=patterns_ch_ref(:,j);

seq_ref_current=patterns_seq_ref(:,j);

ind_occ_ref_current=patterns_ind_occ_ref(j);

%search potential patterns of size s+1 among patterns of size s: the first channel of the pattern

%of reference must be the last one of the pattern of size s AND the duration of the potential new

%pattern must not be > d2 AND the pattern of size s+1 must not already contain d3 times the same

%channel (equivalent to the pattern of size s must not already contain d3-1 times the last channel

%of the pattern of reference)

ind_sub_new=find(ch_ref_current(1)==patterns_ch_current(end,:)&patterns_seq_current(end,:)+...

seq_ref_current(end)<=d2&sum(patterns_ch_current==ch_ref_current(2))<d3);

%occurences times of the current pattern of reference

times_current=patterns_time_ref(ind_occ_ref_current:ind_occ_ref_current+...

patterns_histo_ref(j)-1)-uint32(seq_ref_current(end));

%Optimization: the following instructions allows managing all patterns of size

%s+1 simultaneously

%get the occurences times of patterns of size s previously selected

temp_ind=find(ismember(patterns_ind_pat_current,ind_sub_new));

patterns_time_current_temp=patterns_time_current(temp_ind);

%which pattern do these occurence times refer to

patterns_ind_pat_current_temp=patterns_ind_pat_current(temp_ind);

%select occurences which match those of the pattern of

%reference, i.e. which will be the occurences of the patterns of

%size s+1

ind_intersect=ismember(patterns_time_current_temp,times_current);

%compute the number of occurences for all potential patterns of

%size s+1

histo_new=accumarray(patterns_ind_pat_current_temp',ind_intersect);

histo_new=histo_new(ind_sub_new);

%Apply rule 1

ind=histo_new>d1;

histo_new=histo_new(ind);

ind_sub_new=ind_sub_new(ind);

if ~isempty(ind_sub_new),

%number of patterns of size s+1

nb=length(ind_sub_new);

%store the information for the patterns of size s+1

patterns_ch_seq_new(:,count_patterns:count_patterns+nb-1)=[patterns_ch_current(:,ind_sub_new);...

ch_ref_current(2)*ones(1,nb,'uint8')];

patterns_seq_new(:,count_patterns:count_patterns+nb-1)=[patterns_seq_current(:,ind_sub_new);...

seq_ref_current(end)+patterns_seq_current(end,ind_sub_new)];

patterns_histo_new(count_patterns:count_patterns+nb-1)=histo_new';

patterns_ind_sub_new(count_patterns:count_patterns+nb-1)=ind_sub_new;

%check if storing the new patterns will overflow the current

%matrices. If so, enlarge them.

if count_tot_patterns+sum(histo_new)-1>taille_time,

patterns_time_new=[patterns_time_new zeros(1,taille2,'uint32')];

patterns_time_ind_new=[patterns_time_ind_new zeros(1,taille2,'uint32')];

taille_time=taille_time+taille2;

end;

patterns_time_new(count_tot_patterns:count_tot_patterns+sum(histo_new)-1)=...

patterns_time_current_temp(ind_intersect&ismember(patterns_ind_pat_current_temp,...

ind_sub_new))+uint32(seq_ref_current(end));

%(only works if ind_sub_new is in ascending order)

patterns_time_ind_new(count_tot_patterns:count_tot_patterns+sum(histo_new)-1)=...

count_patterns-1+cumsum([1 diff(patterns_ind_pat_current_temp(ind_intersect&...

ismember(patterns_ind_pat_current_temp,ind_sub_new)))>0]);

count_tot_patterns=count_tot_patterns+sum(histo_new);

count_patterns=count_patterns+nb;

end;

end;

%if some patterns of size s+1 were found

if count_patterns-1>0,

%archive results

clear patterns_time_currentpatterns_ind_pat_current

patterns_ch_archive{s}=patterns_ch_seq_new(:,1:count_patterns-1);

patterns_seq_archive{s}=patterns_seq_new(:,1:count_patterns-1);

patterns_histo_archive{s}=patterns_histo_new(1:count_patterns-1);

patterns_ind_sub_archive{s}=patterns_ind_sub_new(1:count_patterns-1);

fprintf([num2str(s) ' : ' num2str(count_patterns-1) '\n'])

clear patterns_ch_seq_newpatterns_seq_newpatterns_histo_newpatterns_ind_sub_new

%the information for the patterns of size s allowing to find patterns of size s+1 is replaced by

%that for the patterns of size s+1

patterns_ch_current=patterns_ch_archive{s};

patterns_seq_current=patterns_seq_archive{s};

patterns_time_current=patterns_time_new(1:count_tot_patterns-1);

clear patterns_time_new

patterns_ind_pat_current=patterns_time_ind_new(1:count_tot_patterns-1);

clear patterns_time_ind_new

[a,b]=ismember(patterns_time_current,only_times_binned(:,1));

clear a;

if option>0

patterns_special_value{s}=abs(accumarray(patterns_ind_pat_current',...

exp(-2*i*pi*only_times_binned(b,2)*option*binw)));

end;

clear bpatterns_time_ind_new;

else

%no patterns found, stop the loop

non_stop=0;

end;

end;

MI maximization by Backward Method

function [MI_est,remain_col,history_MI_est]=MI_maximization2(distrib,base)

%maximum of mutual information available in a matrix by selecting the

%columns from a backward method

%

%USE :

%[MI_est,remain_col,history_MI_est]=MI_maximization2(distrib,base)

%

%INPUT : distrib - joint distribution between stimuli in rows and observations in

% columns

% base - log base

%

%OUTPUT : MI_est - final MI estimate

% remain_col - remaining columns

% history_MI_est - history of optimization: MI and variable removed

% at each step

%

%Note: distrib is supposed to have no null value

[ncellx,ncelly]=size(distrib);

history_MI_est=zeros(2,ncelly);

if ncelly==0||ncellx==0,

error('Empty distribution');

end;

%margins

margin_y=sum(distrib,1);

margin_x=sum(distrib,2);

count=sum(margin_x);

%estimate mutual information for the initial matrix

norm_x=repmat(margin_x,1,ncelly);

norm_y=repmat(margin_y,ncellx,1);

norm_xy=norm_x.*norm_y;

MI_est=distrib(:).*log(distrib(:)./norm_xy(:))/count;

MI_est(isnan(MI_est))=0;

MI_est=sum(MI_est)+log(count);

MI_est=MI_est/log(base);

%initialize the history of MI maximization

history_MI_est(:,1)=[MI_est;0];

%remaining columns

ind_rest=1:ncelly;

%Optimization: the removing of a column affects margin_x more than margin_y whose

%only one value has to be removed. We then decompose I(X,Y)=H(X)+H(Y)-H(X,Y) and

%compute everything related to margin_y first.

%p(x,y).log(p(x,y))

mat_ref=distrib.*log(distrib);

mat_ref(isnan(mat_ref))=0;

%sum only on x of p(x,y).log(p(x,y))

%equivalent to -H(X,Y) before integration over y

mat_ref_y=sum(mat_ref,1);

%p(y).log(p(y)) for each y

%equivalent to -H(Y) before integration over y

ent_margin_y=margin_y.*log(margin_y);

%H(Y)-H(X,Y) before integration over y

mat_ref_y=mat_ref_y-ent_margin_y;

MI_est=-sum(mat_ref_y)/sum(margin_y)/log(base)-margin_x'/sum(margin_x)*log(margin_x/sum(margin_x))/log(base);

%initialize the history of MI maximization

history_MI_est(:,1)=[MI_est;0];

stop=0;

counter=1;

%while MI continues to increase

while ~stop,

clear MI_i

counter=counter+1;

%compute H(Y)-H(X,Y) for each column removed (sum on y minus what was

%in each y value)

normalization=sum(margin_y(ind_rest))-margin_y(ind_rest);

mat_ref_y_sum=sum(mat_ref_y(ind_rest))-mat_ref_y(ind_rest);

mat_ref_y_sum=mat_ref_y_sum./normalization;

%compute the MI associated to the removal of each column