%  MATLAB tattaiRin.m . Enlarged version of the model of  Thattai- Van 
% Oudenaarden (2001)
% THIS PROGRAM calculates INTRINSIC NOISE for program TATTAIR.m 
%  SINGLE GENE MODEL
%       There are 4 reactions (see text how this arises): 
%   r=mRNA, p=protein, in reactions =
%      stands for right arrow
% 1st reaction is production reaction of mRNA: the effectice rate constant
%  will depend on polymerase R
% reaction 0=r(c1)
% reaction r=0 (c2)
% reaction r=r+p (c3)
% reaction p=0 (c4)
% parameters of the reaction system, nens=nb of cells, timelimit= length
%    of running time, npol = size of loop on number of polymerases,itime 
%  is length of file which gives time 
% dependence of proteins or averages
      itime=1000;
      nens=200;
      npol=32;
      timelimit=40000;
%  parameters:reaction rates
%      b=20;
       b=2;
      c1=0.01;
      c1=2.*c1;
% multiplication of c1 by 2 implements the choice of average behavior for 
% the random variable R (see text)
      c2=log(2)/120;
      c3=b*c2;
      c4=log(2)/3600;
      avp=0;
      avp2=0;
      tim=zeros(itime-1,1);
      npt=zeros(itime-1,1);
      npt2=zeros(itime-1,1);
      np1=zeros(itime-1,1);
      np12=zeros(itime-1,1);
      sint=zeros(itime-1,1);
% parameters for gaussian in R polymerase
      Rav=30;
      mi=30;
      spo=10;
      int1=1/sqrt(2*pi);
      int1=int1/spo;
%  LOOP OVER POLYMERASE
       for ilo=1:npol
      Rp=4+(ilo-1)*2;
%  LOOP OVER CELLS
        for  is=1:nens
%  intialization
      r=1;
      p=0;
      time=0;
      izt=1;
     while time < timelimit
      sum(1)=0;
      sum(2)=0;
      sum(3)=0;
      sum(4)=0;
      sum(5)=0;
      a1=c1.*(Rp./(mi+Rp));
      a2=c2.*r;
      a3=c3.*r;
      a4=c4.*p;
      a0=a1+a2+a3+a4;
      r1=rand(1,1);
      r2=rand(1,1);
      tau=-(1./a0)*log(r1);
      yr2=r2*a0;
      sum(2)=sum(1)+a1;
      sum(3)=sum(2)+a2;
      sum(4)=sum(3)+a3;
      sum(5)=sum(4)+a4;
      for k=2:5
      if (sum(k) >= yr2) & (sum(k-1) < yr2) 
          mu=k-1;
      end 
        end
      if (mu == 1) 
      r=r+1;
      end
      if (mu == 2)
      r=r-1;
       end  
      if (mu == 3) 
      p=p+1;
       end  
      if (mu == 4) 
          p=p-1;
      end
       ztime=timelimit./itime;
       if time > izt*ztime
           time;
           izt;
           izt*ztime;
           np1(izt)=np1(izt)+p;
           np12(izt)=np12(izt)+p^2;
              
           izt=izt+1;
               end 
       time=time + tau;
          end
         end
       imet=int1*exp(-(Rp-Rav)^2./(2*spo^2));
       np1=np1./nens;
       np12=np12./nens;
       coeff=(np12-np1.^2);
       sint=sint+imet.*coeff;
         end
       sint=sqrt(sint);
%  time behavior of proteins averaged over cells
       for jt=1:itime-1
       tim(jt)=ztime*jt;
        end


       tim1=tim([1:itime-10],[1]);
       sint1=sint([1:itime-10],[1]);

       plot(tim1,sint1)
       title('intrinsic noise, averaged over cells, b=2')
       xlabel('time in seconds')