% Order: mn, mc, pc, pn
clear AllTs AllNs
N = [0 0 0 10];
Tmax = 800;
index=1;
V=1;
AllNs(index,:)=N;
AllTs(index)=0;
deltam = 1;
deltap = 0.5;
gammam = 1;
gammap = 1;
alpha = 10;
beta = 1;
r = 2;
a = (alpha*beta)/(deltam*deltap);
k = 0.1;
Results = [1 0 0 0; -1 1 0 0;...
 0 -1 0 0; 0 0 1 0; 0 0 -1 1; 0 0 0 -1];
time=0;
pLer=[];
while time < Tmax
    % Calculate rates
    mn = N(1);
    mc = N(2);
    pc = N(3);
    pn = N(4);
    p0 = 0;
    if ((time > Tmax/4 && time < Tmax/2) || (time > 3*Tmax/4))
        deltap=2;
    else
        deltap = 0.5;
    end
    for ell=0:min(pn,r)
        p0 = p0+(1/k)^ell*factorial(pn)/factorial(pn-ell)*factorial(r)/factorial(r-ell)*1/factorial(ell);
    end
    try
        pLTop = (1/k)^r*factorial(pn)/factorial(pn-r);
    catch
        pLTop = 0;
    end
    pL = pLTop/p0+1e-2;
    pL2 = (pn/(k+pn))^r+1e-2;
    pLer=[pLer;(pL-pL2)/pL];
    Rates = [alpha*pL gammam*mn deltam*mc beta*mc gammap*pc deltap*pn];
    NextTimes = -log(rand(1,length(Rates)))./Rates;
    [tmin,kmin]=min(NextTimes);
    time = time+tmin;
    if (time < Tmax)
        N=N+Results(kmin,:);
        index=index+1;
        AllTs(index)=time;
        AllNs(index,:)=N;
    end
end