% Single trial, second version (for figure 4)
%
% Jorge F. Mejias, 2025
%

function [rate,ratesub]=trial2(par,optopulse,Nareas,Gw)


%we set up the variables for this trial:
bringparam(par);
%cortex:
totalinput=zeros(4,round(triallength/dt),Nareas);
rate=zeros(4,round(triallength/dt),Nareas);transfer=zeros(4,Nareas);
xi=normrnd(0,1,4,round(triallength/dt),Nareas); %noise for re2,ri2,re5,ri5 and both areas
input=zeros(4,Nareas);
%subcortex:
totalinputsub=zeros(7,round(triallength/dt),1); %In order: PC, CN, ZI, VL, Pom, PN, and GrC
ratesub=zeros(7,round(triallength/dt),1);transfersub=zeros(7,1);
xisub=normrnd(0,1,7,round(triallength/dt),1); %noise for subcortex
inputsub=zeros(7,1);Gainsub=[10,1,1,5,0.2,1,1]';
%low-pass cerebellar signal:
xilp1=normrnd(0,1,1,round(triallength/dt),1);
lowpassPC=zeros(1,round(triallength/dt),1);
%low-pass prefrontal signal:
xilp2=normrnd(0,1,1,round(triallength/dt),1);
lowpassPFC=zeros(1,round(triallength/dt),1);
%low-pass external signal (to PC and S1):
xilp3=normrnd(0,1,1,round(triallength/dt),1);
lowpassS1=zeros(1,round(triallength/dt),1);


%first iteration:
rate(:,1,:)=5*(1+tanh(2.*xi(:,1,:))); %between 0 and 10 spikes/s
ratesub(:,1,1)=5*(1+tanh(2.*xisub(:,1,1))); %between 0 and 10 spikes/s
%Now we start the real simulation:
i=2;
for time=2*dt:dt:triallength
	

  %external optostimulation to PC:
  optoPC=0;optoS1=0;optoM1=0;
  stamp1=2;stamp2=2.1;
  %stamp1=0;stamp2=triallength;
  if time>=stamp1 && time<stamp2
      optoPC=optopulse(1);
      optoS1=optopulse(2);
      optoM1=optopulse(3);
  end



  %we compute the input to subcortical areas (in order: PC, CN, ZI, VL, Pom, PN, and GrC):
  inputsub(1,1)=0.1+Jpc*ratesub(7,i-1,1)+optoPC+lowpassPC(1,i-1,1)+lowpassS1(1,i-1,1); %input to PC (Purkinje cells)
  inputsub(2,1)=21-ratesub(1,i-1,1); %input to CN (cerebellar nucleus)
  inputsub(3,1)=ratesub(2,i-1,1)-12.; %input to ZI (zona incerta)
  inputsub(4,1)=ratesub(2,i-1,1)-3.*ratesub(3,i-1,1); %input to VL (ventrolateral thalamus)
  inputsub(5,1)=30+0.2*ratesub(2,i-1,1)-0.5*ratesub(3,i-1,1); %input to Pom (posteromedial complex)
  inputsub(6,1)=Jpn*rate(3,i-1,1)+Jpn*rate(3,i-1,2); %input to PN (pontine nucleus) from S1 and M1
  inputsub(7,1)=Jgrc*ratesub(6,i-1,1); %input to GrC (granule cells)
  totalinputsub(:,i-1,1)=inputsub(:,1);


  %we compute the input to every cortical area:
  %input from thalamus (Pom) to cortical area S1:
  Is1e2=Jpom1*ratesub(5,i-1,1)+optoS1;Is1i2=Jpom2*ratesub(5,i-1,1);
  Is1e5=Jpom3*ratesub(5,i-1,1);Is1i5=Jpom4*ratesub(5,i-1,1);
  %input from thalamus (Pom, VL) to cortical area M1:
  Im1e2=Jpome2*ratesub(5,i-1,1)+optoM1;Im1i2=Jpomi2*ratesub(5,i-1,1);
  Im1e5=Jpome5*ratesub(5,i-1,1)+Jvle5*ratesub(4,i-1,1);Im1i5=Jpomi5*ratesub(5,i-1,1);
  Iext=[Is1e2 Im1e2;Is1i2 Im1i2;Is1e5 Im1e5;Is1i5 Im1i5];

  for k=1:Nareas
    %total input to the k-th area:
    input(:,k)=inputbg+Iext(:,k);
    input(3,2)=input(3,2)+lowpassPFC(1,i-1,1); %transient input from PFC to M1 L5/6E, for movement preparation
    input(1,1)=input(1,1)+lowpassS1(1,i-1,1); %external input to S1 (PC also receives it)
    totalinput(:,i-1,k)=input(:,k)+J*rate(:,i-1,k);
  end
      
  %interareal projections:
  % M1 to S1:
  totalinput(1,i-1,1)=totalinput(1,i-1,1)+Gw*W(1,2,1,2)*rate(3,i-1,2);
  totalinput(1,i-1,1)=totalinput(1,i-1,1)+Gw*W(1,2,1,1)*rate(1,i-1,2);
  totalinput(3,i-1,1)=totalinput(3,i-1,1)+Gw*W(1,2,2,2)*rate(3,i-1,2);
  totalinput(2,i-1,1)=totalinput(2,i-1,1)+Gw*W(1,2,3,2)*rate(3,i-1,2);
  totalinput(4,i-1,1)=totalinput(4,i-1,1)+Gw*W(1,2,4,2)*rate(3,i-1,2);
  % S1 to M1:
  totalinput(1,i-1,2)=totalinput(1,i-1,2)+Gw*W(2,1,1,1)*rate(1,i-1,1);
  totalinput(3,i-1,2)=totalinput(3,i-1,2)+Gw*W(2,1,2,1)*rate(1,i-1,1);

  
  for k=1:Nareas
    %input after transfer functions:
    for j=1:4 %for cortex:
      transfer(j,k)=totalinput(j,i-1,k)/(1-exp(-totalinput(j,i-1,k)));
    end
    for j=1:7 %for subcortex:
        transfersub(j,1)=Gainsub(j,1).*totalinputsub(j,i-1,1)/(1-exp(-15.*totalinputsub(j,i-1,1))); %beta=15 is good aprox
    end
  end

  
  %we update the firing rates of all areas:
  for k=1:Nareas
    rate(:,i,k)=rate(:,i-1,k)+tstep.*...
    (-rate(:,i-1,k)+transfer(:,k))+tstep2.*xi(:,i-1,k);
    ratesub(:,i,1)=ratesub(:,i-1,1)+tstep(1).*...
    (-ratesub(:,i-1,1)+transfersub(:,1))+0.5*tstep2(1).*xisub(:,i-1,1);
    lowpassPC(1,i,1)=lowpassPC(1,i-1,1)+(dt/taulpPC)*(-lowpassPC(1,i-1,1))+sqrt(dt/taulpPC)*siglpPC*xilp1(1,i-1,1);
    lowpassPFC(1,i,1)=lowpassPFC(1,i-1,1)+(dt/taulpPFC)*(-lowpassPFC(1,i-1,1))+sqrt(dt/taulpPFC)*siglpPFC*xilp2(1,i-1,1);
    lowpassS1(1,i,1)=lowpassS1(1,i-1,1)+(dt/taulpS1)*(-lowpassS1(1,i-1,1))+sqrt(dt/taulpS1)*siglpS1*xilp3(1,i-1,1);
  end
  
  
  %index iteration
  i=i+1;
end