function [ mout,pout,rout ] = GillespieTiming(N,parms,totalreps,mstart, pstart, rstart, tau,Output_Times) 
%DELAYGILLESPIE Runs Gillespie algorithm with delay processes. Returns a
% value of time (dist)
%unpack parms
P0=parms.p0;
X1=parms.x1;
X2=parms.x2;
P1=parms.p1;
P2=parms.p2;
P3=parms.p3;
NP=parms.n;
N1=parms.n1;
N2=parms.n2;
N3=parms.n3;
MUP=parms.mup;
ALPHAP=parms.alphap;
ALPHAR=parms.alphar;
ALPHAM=parms.alpham;
MURR=parms.murr;

mout = zeros(totalreps,length(Output_Times));
pout = zeros(totalreps,length(Output_Times));
rout = zeros(totalreps,length(Output_Times));


for reps=1:totalreps %parfor
    
j_t_next = 1;

    
iter = 1;
rlist = [];

% sets initial values
m=mstart;
mn=0;
p=pstart;
pn=0;
r=rstart;
rn=0;
t=0;tn=0;
rlist=[];

%When rlist is empty, there is a different process:
%We follow this process until rlist is no longer empty
a1 = N*(X1 + (X2-X1)/(1+((r/N)/P1)^N1))*(m/N);
a2 = MUP*p; 
a3 = N*(ALPHAP*(m/N))/(1 + ((r/N)/P2)^N2);
a4 = N*ALPHAR/(1 + ((p/N)/P3)^N3);
a5 = MURR*(r);
a6 = N*ALPHAM/(1 + ((p/N)/P0)^NP);
while t<Output_Times(end) % runs until p hits a threshold            
            a0=a1+a2+a3+a4+a5+a6;
            r1=rand(1);
            r2=rand(1);
            dt=(1/a0)*log(1/r1);
            if numel(rlist)>0 && t<=rlist(1) && rlist(1)<=(t+dt)
            %if le(t(iter),rlist(1))*le(rlist(1),(t(iter)+dt))
                mn= m+1;
                pn = p;
                rn = r;
                tn = rlist(1);
                rlist(1)=[];
                a1=(X1 + (X2-X1)/(1+((rn/N)/P1)^N1))*(mn);
                a3=N*(ALPHAP*(mn/N))/(1 + ((rn/N)/P2)^N2);
                a5 = MURR*(rn);
            else
                if r2*a0<=a1
                    mn= m-1;
                    pn = p;
                    rn = r;
                    a1=(X1 + (X2-X1)/(1+((rn/N)/P1)^N1))*(mn);
                    a3=N*(ALPHAP*(mn/N))/(1 + ((rn/N)/P2)^N2);
                    a5 = MURR*(rn);
                %elseif le(a1,r2*a0)*le(r2*a0,(a1+a2))
                elseif a1<=r2*a0 && r2*a0<=(a1+a2)
                    mn = m;
                    pn = p-1;
                    rn = r;
                    a2=MUP*pn;
                    a4 = N*ALPHAR/(1 + ((pn/N)/P3)^N3);
                    a6=N*ALPHAM/(1 + ((pn/N)/P0)^NP);
                elseif (a1+a2)<=r2*a0 && r2*a0<=(a1+a2+a3)
                    mn = m;
                    pn = p+1;
                    rn = r;
                    a2=MUP*pn;
                    a4 = N*ALPHAR/(1 + ((pn/N)/P3)^N3);
                    a6=  N*ALPHAM/(1 + ((pn/N)/P0)^NP);
                elseif (a1+a2+a3)<=r2*a0&& r2*a0<=(a1+a2+a3+a4)
                    mn = m;
                    pn = p;
                    rn = r+1;
                    a1=(X1 + (X2-X1)/(1+((rn/N)/P1)^N1))*(mn); 
                    a3=N*(ALPHAP*(mn/N))/(1 + ((rn/N)/P2)^N2);
                    a5 = MURR*(rn);
                elseif (a1+a2+a3+a4)<=r2*a0 && r2*a0<=(a1+a2+a3+a4+a5)
                    mn = m;
                    pn = p;
                    rn = r-1;
                    a1=(X1 + (X2-X1)/(1+((rn/N)/P1)^N1))*(mn); 
                    a3=N*(ALPHAP*(mn/N))/(1 + ((rn/N)/P2)^N2);
                    a5 = MURR*(rn);
                else
                    rlist = [rlist (t+tau)];
                    mn = m;
                    pn = p;
                    rn = r;
                end
                tn = t+dt;
            end
        iter = iter + 1;
        t=tn;
        m =mn;
        p =pn;
        r =rn;
        while j_t_next<=length(Output_Times)&&t>Output_Times(j_t_next)
            mout(reps,j_t_next) = m;
            pout(reps,j_t_next) = p;
            rout(reps,j_t_next) = r;
            j_t_next=j_t_next+1;
        end
end 
   
end                   
                                     
end

