function Figure5_1b

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

k1=0.18;
k2=0.00025;
k3=37.5;
k4=2200;

X=0;
a0=k1*X*(X-1)+k2*X*(X-1)*(X-2)+k3*X+k4;

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

while (time<100)
   timeSSAprev=timeSSA; 
   timefinSSA=500;
   kk=kk+1;
   while (timeSSA<timeSSAprev+timefinSSA)
      rr=rand(2,1);
      timeSSA=timeSSA+1;
      a0=k1*X*(X-1)+k2*X*(X-1)*(X-2)+k3*X+k4;
      tau=(1/a0)*log(1/rr(1));
      ss=k1*X*(X-1);
      if (ss>rr(2)*a0)
          X=X+1;
      else
          ss=ss+k2*X*(X-1)*(X-2);
          if (ss>rr(2)*a0)
             X=X-1;
          else
             ss=ss+k3*X;     
             if (ss>rr(2)*a0)
                 X=X-1;
             else
                 ss=ss+k4;
                 if (ss>rr(2)*a0)
                    X=X+1;
                 else
                    STOP
                 end;
             end;    
          end;       
      end;           
      time=time+tau;
   end;     
XX(kk)=X;
tt(kk)=time;
end;
[t,z] = ode45(@myode,[0 20],[0]);

figure(1);
set(gca,'Fontsize',20);
plot(tt,XX,'b','Linewidth',1);
hold on;
plot(t,z,'r','Linewidth',4);
line([20 100] ,[100 100],'Color','r','Linewidth',4);
xlabel('time [min]');
ylabel('number of molecules');
axis([0 100 0 550]);

function dydt = myode(t,z)
k1=0.18;
k2=0.00025;
k3=37.5;
k4=2200;
dydt = [k1*z(1)*z(1)-k2*z(1)*z(1)*z(1)-k3*z(1)+k4];