%
% Large-scale model of the sensory trigemino-thalamo-cerebelo-cortical
% circuit, from the paper Romano et al. Nature Communications (2025)
%
%
% Model and code by Jorge Mejias, 2025
%
% -------------------


format short;clear all;
close all;clc;rng(938199);


%PFCsig reflects prefrontal signals during movement preparation, leading the PC-M1 correlation (so lag(M1,PC)<0).
%So, for premovement condition use PFCsig=100. For movement condition, use PFCsig=0;
optopulse=[0 0 0]; %pulse strengths for PC, VL/Pom and PN, respectively. Use 300 for VL/Pom/PN, and 0.6 for PC
PCsig=1*0.5;sigext=1*2.;PFCsig=0*100; %PCsig=[0,0.5], sigext=[0,2], PFCsig=[0,100].
condition='Pre-movement';if PFCsig<1;condition='Movement';end


%First, let's plot an example of cross-correlation plots:
[par,rate,ratesub,S1,M1,corrS,lagsS,corrM,lagsM]=conditionsim(optopulse,PCsig,PFCsig,sigext);
ratePC=ratesub(1,:,1);corrlim=100;multifigures1;
lagsplot=lagsM;corrplot=corrM;


%now let's do a bit more statistics:
Ntrials=10;
lagdistS=zeros(1,Ntrials);corrdistS=lagdistS;
lagdistM=zeros(1,Ntrials);corrdistM=lagdistM;
pearsonS=zeros(1,Ntrials);pearsonM=pearsonS;

Tlength=round(par.triallength/par.dt);
fulltraces=zeros(Ntrials,Tlength,3); %trials for PC, S1, M1

for i=1:Ntrials
   %First, simulate:
   [par,rate,ratesub,S1,M1,corrS,lagsS,corrM,lagsM]=conditionsim(optopulse,PCsig,PFCsig,sigext);
   ratePC=ratesub(1,:,1);corrlim=100;

   %we save the traces for each trial:
   fulltraces(i,:,1)=ratePC;
   fulltraces(i,:,2)=S1;
   fulltraces(i,:,3)=M1;

   %Let's compute the peak cross-correlation and its time lag for PC-M1 correlation
   Lags=lagsM*par.dt*1000;
   [~,b]=find(abs(Lags)<corrlim); %'b' has the indices for lags between -corrlim and corrlim.
   [c,d]=max(corrM(b));corrdistM(1,i)=c;lagdistM(1,i)=-corrlim+d*par.dt*1000;

   %Now let's compute the peak cross-correlation and its time lag for PC-S1 correlation
   Lags=lagsS*par.dt*1000;
   [~,b]=find(abs(Lags)<corrlim); %'b' has the indices for lags between -corrlim and corrlim.
   [c,d]=max(corrS(b));corrdistS(1,i)=c;lagdistS(1,i)=-corrlim+d*par.dt*1000;

   %Now, let's compute the Pearson correlation coefficient between PC and M1:
   sampling=800;dPC=downsample(ratePC,sampling);dM1=downsample(M1,sampling);
   [rho,pval]=corrcoef(dM1,dPC);pearsonM(1,i)=rho(1,2);

   %And between PC and S1:
   sampling=800;dPC=downsample(ratePC,sampling);dS1=downsample(S1,sampling);
   [rho,pval]=corrcoef(dS1,dPC);pearsonS(1,i)=rho(1,2);

end

%%%%%cross-correlations:
disp('Cross-correlations (lag and peak, first PC-S1 then PC-M1):')
mean(lagdistS)
mean(corrdistS)

mean(lagdistM)
mean(corrdistM)

% Is the peak correlation between PC-S1 significantly larger that the one 
% between PC-M1? If pSM<0.05, then yes.
disp('Significant?');[~,pSM]=ttest(corrdistS,corrdistM)

comparisons=[lagdistS' corrdistS' pearsonS' lagdistM' corrdistM' pearsonM'];

%%%%Pearson correlations:
disp('Pearson correlations (first PC-S1, then PC-M1):')
mean(pearsonS)
mean(pearsonM)
disp('Significant?');[~,pSM]=ttest(pearsonS,pearsonM)


save traces.mat fulltraces comparisons pSM optopulse PFCsig par lagsplot corrplot