% Loop of four reactions
Nspec = 4;
DeltaNs = zeros(2*Nspec,Nspec);
for iN=1:Nspec
    % Forward reactions first
    DeltaNs(iN,iN)=-1;
    DeltaNs(iN,modn(iN+1,Nspec))=1;
    % Backward reactions
    DeltaNs(Nspec+iN,iN)=1;
    DeltaNs(Nspec+iN,modn(iN+1,Nspec))=-1;
end
rng(0);
ForwardRates = rand(1,4);
BackwardRates = rand(1,4);
% Detailed balance
ForwardRates(1) = prod(BackwardRates)/prod(ForwardRates(2:end));
N0 = 1000;
Ns = [N0 0 0 0];
allNs = Ns;
allts = 0;
Tmax = 100;
t=0;
while (t < Tmax)
    Rates = [ForwardRates BackwardRates].*[Ns Ns(2:end) Ns(1)];
    Times = -log(rand(1,2*Nspec))./Rates;
    [tmin,kmin]=min(Times);
    t=t+tmin;
    if (t < Tmax)
        Ns = Ns + DeltaNs(kmin,:);
        allNs = [allNs; Ns];
        allts = [allts; t];
    end
end
% Plot the expected mean 
MeanMat = zeros(4);
MeanMat(1,:)=1;
for iN=1:Nspec-1
    MeanMat(iN+1,iN)=ForwardRates(iN);
    MeanMat(iN+1,iN+1)=-BackwardRates(iN);
end
Xbar = MeanMat \ [N0;zeros(Nspec-1,1)];
plot(allts,allNs')
set(gca,'ColorOrderIndex',1)
hold on
plot(xlim,[Xbar Xbar],':')
% Estimate fluxes
startN = find(allts > 10);
startN = startN(1);
Fluxes=[MeanMat(2:end,:); -BackwardRates(4) 0 0 ForwardRates(4)]*mean(allNs(startN:end,:))'

function x = modn(a,b)
    x = mod(a,b);
    if (x==0)
        x=b;
    end
end