%
% Simulations of the system when we transiently stimulate different
% regions (such as PN, or thalamus, or cortex.
%
% Jorge Mejias, 2025
%


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


%Simulate the full system, with all sigmas at zero:
%PFCsig should reflect the prefrontal signals during movement preparation,
%which should lead to cortex leading the PC-M1 correlation (so lag(M1,PC)<0).
%So, for premovement condition use PFCsig=100. For movement condition, use PFCsig=0;
optopulse=[1*0.6 0*60 0*60]; %Pulses for PC, S1 and M1. Use 0.6 for PC, and 60 for S1/M1.
PCsig=0*0.5;sigext=0*2;PFCsig=0*100.; %all sigmas set to zero to see the effect more clearly.
Nareas=2;Gw=1.;Tmax=5.;par=parameters(Tmax,PCsig,PFCsig,sigext);
condition='Pre-movement';if PFCsig<1;condition='Movement';end
Ntrials=10;S1=zeros(Ntrials,round(par.triallength/par.dt));M1=S1;PC=S1;


%simulate all the trials:
for i=1:Ntrials
    [rate,ratesub]=trial2(par,optopulse,Nareas,Gw);ratePC=ratesub(1,:,1);eta=0.2;
    PC(i,:)=ratesub(1,:,1);
    S1(i,:)=eta*rate(1,:,1)+(1-eta)*rate(3,:,1);
    M1(i,:)=eta*rate(1,:,2)+(1-eta)*rate(3,:,2);
end
time=par.dt:par.dt:Tmax;


%We plot the system's response to a pulse (for clarity, set all three sigmas above to zero):
cmap=colororder();Ntrialsplotted=Ntrials;
for kk=1:Ntrialsplotted
    plot(time,5+PC(kk,:),'LineWidth',2,'Color',cmap(1,:));hold on;plot(time,15+S1(kk,:),'-','LineWidth',2,'Color',cmap(2,:));hold on;plot(time,15+M1(kk,:),'-','LineWidth',2,'Color',cmap(3,:));hold on; 
end

set(gca,'FontSize',12,'LineWidth',2,'TickLength',[0.02 0.02]);
set(gca,'box','off');xlabel('Time (s)');ylabel('Firing rate (sp/s)');set(gca,'Layer','top');
legend({'PC','S1','M1'});xlim([1.9 2.3]);ylim([4 40]);

%save stimPC.mat Ntrials time PC S1 M1