Noei_PNAS_2022/Main/corr_sim.m

%% routine for the simulation to show NMF vs correlation and generating figure 1


%% adding necessary directories and loading data
clear all 
clc 
close all


filepath = fileparts(pwd);
addpath(strcat(filepath,filesep,'Toolboxes',filesep,'PopulationSpikeTrainFactorizationV1.0',filesep));
addpath(strcat(filepath,filesep,'Aux_functions',filesep));
addpath(strcat(filepath,filesep,'Toolboxes',filesep,'BCT',filesep,'2019_03_03_BCT',filesep));




%% Simulations to show examples of correlation vs NMF
n_neurons = 10;
n_bins = 60*10;
base_line_rate = 1;
n_permutes = 100;

% case 1
correlated_neurons = 2:7;
correlated_time = 1:300;
SNR = 1.5;
rates = 0.1*randn(n_neurons,n_bins) + base_line_rate;
rates(correlated_neurons,correlated_time) = 0.1*rand(length(correlated_neurons),length(correlated_time)) + SNR*base_line_rate;
rates(9:10,end/2:end) = rates(9:10,end/2:end) + 0.5;

C_real = corr(rates');
corr_sh = zeros(n_neurons,n_neurons,n_permutes);
for permute=1:n_permutes
    shuffled_spikes = rates;
    for unit=1:n_neurons
        shuffled_spikes(unit,:) = rates(unit,randperm(n_bins));
    end
    corr_sh(:,:,permute) = corr(shuffled_spikes');
end
th_corr = prctile(corr_sh,99,3);
th_corr(th_corr < 0) = 0;
C_th  = C_real.*(C_real > th_corr);
[M,Q] = community_louvain(C_th);

var_explained = zeros(1,n_neurons/2);
for i=1:n_neurons/2
    N_modules = i;
    [~,~,err] = nmf(rates',N_modules);
    var_explained(i) = 1 - (err./norm(rates,'fro'));
end
[idxOfBestPoint] = elbow_finder(var_explained);
[W,H] = nmf(rates',idxOfBestPoint,[],2,500,[],[],5);
H_norm = normalize(H,2,'range');
W_norm = normalize(W,'range');



subplot(3,4,1);imagesc(rates);colormap gray; colorbar;axis square;title('spike rates');
ylabel('Neurons');caxis([1,3])
subplot(3,4,2);imagesc(C_th>0);colormap gray; colorbar;axis square;caxis([0,1]);
plot_comunity(M);title('thresholded correlation');
subplot(3,4,3);imagesc(H_norm');colormap gray;axis square;colorbar;title('spatial modules');
subplot(3,4,4);imagesc(W_norm');colormap gray;axis square;colorbar;title('temporal modules');

% case 2
rates = 0.1*rand(n_neurons,n_bins) + base_line_rate;
correlated_neurons_1 = 2:4;
correlated_neurons_2 = 5:7;
correlated_time_1 = 1:300/4;
correlated_time_2 = correlated_time_1(end):300/2;
correlated_time_3 = correlated_time_2(end):300;
SNR = 2;
rates(correlated_neurons_1,correlated_time_1) = 0.1*randn(length(correlated_neurons_1),length(correlated_time_1)) + SNR*base_line_rate;
rates(correlated_neurons_2,correlated_time_2) = 0.1*randn(length(correlated_neurons_2),length(correlated_time_2)) + SNR*base_line_rate;
rates(correlated_neurons,correlated_time_3) = 0.1*randn(length(correlated_neurons),length(correlated_time_3)) + 1.1*SNR*base_line_rate;
rates(9:10,end/2:end) = rates(9:10,end/2:end) + 0.5;

C_real = corr(rates');
corr_sh = zeros(n_neurons,n_neurons,n_permutes);
for permute=1:n_permutes
    shuffled_spikes = rates;
    for unit=1:n_neurons
        shuffled_spikes(unit,:) = rates(unit,randperm(n_bins));
    end
    corr_sh(:,:,permute) = corr(shuffled_spikes');
end
th_corr = prctile(corr_sh,99,3);
th_corr(th_corr < 0) = 0;
C_th  = C_real.*(C_real > th_corr);
[M,Q] = community_louvain(C_th);

var_explained = zeros(1,n_neurons/2);
for i=1:n_neurons/2
    N_modules = i;
    [~,~,err] = nmf(rates',N_modules);
    var_explained(i) = 1 - (err./norm(rates,'fro'));
end
[idxOfBestPoint] = elbow_finder(var_explained);
[W,H] = nmf(rates',idxOfBestPoint,[],2,500,[],[],5);
H_norm = normalize(H,2,'range');
W_norm = normalize(W,'range');



subplot(3,4,5);imagesc(rates);colormap gray; colorbar;axis square;
ylabel('Neurons');caxis([1,3])
subplot(3,4,6);imagesc(C_th>0);colormap gray; colorbar;axis square;caxis([0,1]);
plot_comunity(M);
subplot(3,4,7);imagesc(H_norm');colormap gray;axis square;colorbar
subplot(3,4,8);imagesc(W_norm');colormap gray;axis square;colorbar


% case 3
correlated_neurons = 2:7;
correlated_time_1 = 1:150;
correlated_time_2 = n_bins-150:n_bins;
base_line_rate = 1;
SNR = 1.5;
rates = 0.1*randn(n_neurons,n_bins) + base_line_rate;
rates(correlated_neurons,correlated_time_1) = 0.1*randn(length(correlated_neurons),length(correlated_time_1)) + SNR*base_line_rate;
rates(correlated_neurons,correlated_time_2) = 0.1*randn(length(correlated_neurons),length(correlated_time_2)) + SNR*base_line_rate;
rates(9:10,end/2:end) = rates(9:10,end/2:end) + 0.5;
C_real = corr(rates');
corr_sh = zeros(n_neurons,n_neurons,n_permutes);
for permute=1:n_permutes
    shuffled_spikes = rates;
    for unit=1:n_neurons
        shuffled_spikes(unit,:) = rates(unit,randperm(n_bins));
    end
    corr_sh(:,:,permute) = corr(shuffled_spikes');
end
th_corr = prctile(corr_sh,99,3);
th_corr(th_corr < 0) = 0;
C_th  = C_real.*(C_real > th_corr);
[M,Q] = community_louvain(C_th);

var_explained = zeros(1,n_neurons/2);
for i=1:n_neurons/2
    N_modules = i;
    [~,~,err] = nmf(rates',N_modules);
    var_explained(i) = 1 - (err./norm(rates,'fro'));
end
[idxOfBestPoint] = elbow_finder(var_explained);
[W,H] = nmf(rates',idxOfBestPoint,[],2,500,[],[],5);
H_norm = normalize(H,2,'range');
W_norm = normalize(W,'range');


subplot(3,4,9);imagesc(rates);colormap gray; colorbar;axis square
ylabel('Neurons');xlabel('Time (ms)');caxis([1,3])
subplot(3,4,10);imagesc(C_th>0);colormap gray; colorbar;axis square;caxis([0,1]);
plot_comunity(M);
subplot(3,4,11);imagesc(H_norm');colormap gray;axis square;colorbar
subplot(3,4,12);imagesc(W_norm');colormap gray;axis square;colorbar

%% Generating the figure relating this concept to cortical Rhythms
% type 1
n_neurons = 10;
n_bins = 60*10;

base_line_rate = 0.1;

rates = 0.001*rand(n_neurons,n_bins) + base_line_rate;
correlated_neurons = 2:7;

correlated_time = 1:300;
SNR = 2;
rates(correlated_neurons,correlated_time) = rand(length(correlated_neurons),length(correlated_time)) + 1.1*SNR*base_line_rate;
rates(9:10,end/2:end) = rates(9:10,end/2:end) + 1;

to_raster = cell(1,n_neurons);
for neuron=1:n_neurons
    to_raster{neuron} = PechePourPoisson(rates(neuron,:),0.1)';
end

t = 0:0.01:60;
LFP = cos(2*pi*t/10) + cos(pi*t/10) + cos(pi*t/100); 
LFP(1:3000) = LFP(1:3000) - 2.5*cos(0.3*pi*t(1:3000)/2);



figure();
subplot(2,1,1);plotSpikeRaster(to_raster,'plottype','vertline');
subplot(2,1,2);plot(t(1:3000),LFP(1:3000));hold all
plot(t(3001:end),LFP(3001:end));


% type 2
n_neurons = 10;
n_bins = 60*10;

base_line_rate = 0.1;

rates = 0.001*rand(n_neurons,n_bins) + base_line_rate;
correlated_neurons = 2:7;
correlated_neurons_1 = 2:4;
correlated_neurons_2 = 5:7;
correlated_time_1 = 1:300/4;
correlated_time_2 = correlated_time_1(end):300/2;
correlated_time_3 = correlated_time_2(end):300;
SNR = 2;
rates(correlated_neurons_1,correlated_time_1) = rand(length(correlated_neurons_1),length(correlated_time_1)) + 0.9*SNR*base_line_rate;
rates(correlated_neurons_2,correlated_time_2) = rand(length(correlated_neurons_2),length(correlated_time_2)) + 0.9*SNR*base_line_rate;
rates(correlated_neurons,correlated_time_3) = rand(length(correlated_neurons),length(correlated_time_3)) + 1.1*SNR*base_line_rate;
rates(9:10,end/2:end) = rates(9:10,end/2:end) + 1;

to_raster = cell(1,n_neurons);
for neuron=1:n_neurons
    to_raster{neuron} = PechePourPoisson(rates(neuron,:),0.1)';
end

t = 0:0.01:60;
LFP = cos(2*pi*t/10) + cos(pi*t/10) + cos(pi*t/100); 
LFP(1:750) = LFP(1:750) - 1.5*cos(8*pi*t(1:750)/10);
LFP(751:1500) = LFP(751:1500) + 2*cos(4*pi*t(751:1500)/10);
LFP(1501:3000) = - 1.5*cos(8*pi*t(1501:3000)/10) + 2*cos(4*pi*t(1501:3000)/10) + LFP(1501:3000);

LFP(751) = mean([LFP(750),LFP(752)]);
LFP(1501) = mean([LFP(1500),LFP(1502)]);


figure();
subplot(2,1,1);plotSpikeRaster(to_raster,'plottype','vertline');
subplot(2,1,2);plot(t(1:750),LFP(1:750));hold all
plot(t(751:1500),LFP(751:1500));plot(t(1501:3000),LFP(1501:3000));
plot(t(3001:end),LFP(3001:end));


% type 3
n_neurons = 10;
n_bins = 60*10;

base_line_rate = 0.1;

rates = 0.001*rand(n_neurons,n_bins) + base_line_rate;
correlated_neurons = 2:7;

correlated_time_1 = 1:300/4;
correlated_time_2 = 3*n_bins/4:n_bins;

SNR = 2;
rates(correlated_neurons,correlated_time_1) = rand(length(correlated_neurons),length(correlated_time_1)) + 0.9*SNR*base_line_rate;
rates(correlated_neurons,correlated_time_2) = rand(length(correlated_neurons),length(correlated_time_2)) + 0.9*SNR*base_line_rate;
rates(9:10,end/2:end) = rates(9:10,end/2:end) + 1;

to_raster = cell(1,n_neurons);
for neuron=1:n_neurons
    to_raster{neuron} = PechePourPoisson(rates(neuron,:),0.1)';
end

t = 0:0.01:60;
LFP = cos(2*pi*t/10) + cos(pi*t/10) + cos(pi*t/100); 
LFP(1:750) = LFP(1:750) - 2*cos(4*pi*t(1:750)/10);
LFP(750:3000) = 2*sin(pi*t(750:3000)/10);
LFP(4500:end) = LFP(4500:end) + 3.5*sin(8*pi*t(4500:end)/10) - 2*cos(4*pi*t(4500:end)/10); 


figure();
subplot(2,1,1);plotSpikeRaster(to_raster,'plottype','vertline');
subplot(2,1,2);plot(t(1:750),LFP(1:750));hold all
plot(t(751:3000),LFP(751:3000));
plot(t(3001:4500),LFP(3001:4500));
plot(t(4500:end),LFP(4500:end));