function Figure5_2

close all;
rand('state',5);

k1=50;
k2=10;
k3=0.00004;
k4=25;

X=10;
Y=10;
a0=k1+k2*X+k3*X*(X-1)*Y+k4;

time=0;
timeSSA=0;
kk=0;

while (time<5000)
   timeSSAprev=timeSSA; 
   timefinSSA=500;
   kk=kk+1;
   while (timeSSA<timeSSAprev+timefinSSA)
      if (X>100) 
          timefinSSA=10;
      end;    
      rr=rand(2,1);
      timeSSA=timeSSA+1;
      a0=k1+k2*X+k3*X*(X-1)*Y+k4;
      tau=(1/a0)*log(1/rr(1));
      ss=k1;
      if (ss>rr(2)*a0)
          X=X+1;
      else
          ss=ss+k2*X;
          if (ss>rr(2)*a0)
             X=X-1;
          else
             ss=ss+k3*X*(X-1)*Y;     
             if (ss>rr(2)*a0)
                 X=X+1;
                 Y=Y-1;
             else
                 ss=ss+k4;
                 if (ss>rr(2)*a0)
                    Y=Y+1;
                 else
                    STOP
                 end;
             end;    
          end;       
      end;           
      time=time+tau;
   end;     
XX(kk)=X;
YY(kk)=Y;
tt(kk)=time;
end;

[t,z] = ode45(@myode,[0 5000],[10; 10]);

figure(1);
set(gca,'Fontsize',20);
h=semilogy([100 100],[20000 20000]);
set(h,'Color','b','Linewidth',4);
hold on;
h=semilogy([100 100],[20000 20000]);
set(h,'Color','r','Linewidth',4);
h=semilogy(tt/60,XX);
set(h,'Color','b','Linewidth',2);
h=semilogy(t/60,z(:,1));
set(h,'Color','r','Linewidth',4);
xlabel('time [min]');
ylabel('number of A molecules');
legend('stochastic','deterministic');
axis([0 82 1 10000]);

figure(2);
set(gca,'Fontsize',20);
h=semilogx([20000 20000],[20000 20000]);
set(h,'Color','b','Linewidth',4);
hold on;
h=semilogx([20000 20000],[20000 20000]);
set(h,'Color','r','Linewidth',4);
xval1=[2:0.15:24.8 25:1:100 110:10:1000 1100:100:6200];
xval2=[2:0.15:24.8 25:1:160];
ynul1=(k2*xval1-k1)./(k3*xval1.*xval1);
ynul2=k4./(k3*xval2.*xval2);
semilogx(XX,YY);
h=semilogx(z(1,1),z(1,2));
set(h,'Color','r','Linewidth',2);
xlabel('number of A molecules');
ylabel('number of B molecules');
h=semilogx(xval1,ynul1);
set(h,'Color','g','Linewidth',4);
h=semilogx(xval2,ynul2);
set(h,'Color','g','Linewidth',4);
h=semilogx([(k4+k1)/k2],[k4*k2*k2/(k3*(k4+k1)*(k4+k1))],'ro');
set(h,'Markersize',4,'Linewidth',6);
axis([0 10000 0 13000]);
h=semilogx(z(:,1),z(:,2));
set(h,'Color','r','Linewidth',2);
legend('stochastic','deterministic');

function dydt = myode(t,z)
k1=50;
k2=10;
k3=0.00004;
k4=25;
dydt =[k1-k2*z(1)+k3*z(1)*z(1)*z(2); k4-k3*z(1)*z(1)*z(2)];

