XcIrna = 1.2;
XcIprot = 1.2;
croRNA = .8;
croPROT=.8;
omegacI = 50;
nucI = 50;
omegacro=50;
nucro=50;
KcI=10;
Kcro=10;
maxn=50000;
x=zeros(1,maxn+1);
y=zeros(1,maxn+1);
z=zeros(1,maxn+1);
w=zeros(1,maxn+1);
v=max(1,8);  %rates
vacc=max(1,8); %rate accumulation
time=zeros(1,maxn+1);

for i=1:50000
    v(1)=omegacI*y(i);
    v(2)=XcIprot*x(i);
    v(3)=nucI*(1-(z(i).^2 / (Kcro.^2 + z(i).^2)));
    v(4)=XcIrna*y(i);
    v(5)=omegacro*w(i);
    v(6)=croPROT*z(i);
    v(7)=nucro*(1-(x(i).^2 / (KcI.^2 + x(i).^2)));
    v(8)=croRNA*w(i);
    R=sum(v);
    y = rand();
    tau = -log(y)/(R);
    time(i+1) = time(i)+tau;

    for m=1:8
        vacc(m)=sum(v(1:m));
    end

    r1=rand();
    compareval=r1*R;
    F=find(vacc>compareval);
    next_rxn=F(1);
    if (next_rxn == 1)
        x(i+1)= x(i) +1;
    elseif (next_rxn == 2)
        x(i+1)= x(i) -1;
    else
        x(i+1)= x(i);
    end

    if (next_rxn == 3)
        y(i+1)=y(i)+1;
    elseif (next_rxn ==4)
        y(i+1)=y(i)-1;
    else
        y(i+1)=y(i);
    end

    if (next_rxn==5)
        z(i+1)=z(i)+1;
    elseif (next_rxn ==6)
        z(i+1)=z(i)-1;
    else
        z(i)=z(i);
        if (next_rxn==7)
            w(i+1)=w(i)+1;
        elseif (next_rxn==8)
            w(i+1)=w(i)-1;
        else
            w(i)=w(i);
        end
    end
end