function [ Dist1,Dist2] = CLEoptTIMINGwhileBINOMIAL( tau,v,Tend,timestep,parms,totalreps )
%EULERFUNCTION Summary of this function goes here
%   Detailed explanation goes here

ttot = (Tend+tau)/timestep;
h = tau/timestep;
mup=parms.mup;
n=parms.n;  
p0=parms.p0;  
a1=parms.a1; 
a2=parms.a2;
p1=parms.p1; 
murr=parms.murr;
alphar=parms.alphar;
alphap=parms.alphap;
alpham=parms.alpham;
p2=parms.p2; 
n2=parms.n2;
n1=parms.n1;
p3=parms.p3; 
n3=parms.n3;
N=parms.N;

vrnd1 = binornd((round(v(4))*ones(1,totalreps)),0.5);
vrnd2 = binornd((round(v(5))*ones(1,totalreps)),0.5);
vrnd3 = binornd((round(v(6))*ones(1,totalreps)),0.5);


parfor j=1:totalreps
    
y = zeros(6,ttot);
vcurr = v;
vcurr(1) = vrnd1(j)/N;
vcurr(2) = vrnd2(j)/N;
vcurr(3) = vrnd3(j)/N;
vcurr(4) = vrnd1(j);
vcurr(5) = vrnd2(j);
vcurr(6) = vrnd3(j);
hist = repmat(vcurr,1,h);
y(:,1:h)=hist;
i = h+1;  

while i<ttot


ylag1 = y(2,(i-(h)) );
ylag2 = y(5,(i-(h)) );
m=y(1,(i-1));
p=y(2,(i-1));
r=y(3,(i-1));

M=y(4,(i-1))/N;
P=y(5,(i-1))/N;
R=y(6,(i-1))/N;
ylag2 = ylag2/N;

Gx = 1/(1+(ylag2/p0)^n);
Grx = 1/(1+(P/p3)^n3);
Fx = 1/(1+(R/p2)^n2);
Sx=(a1)+(a2-a1)/(1+((R)/p1)^n1);


rnds = randn(3,1);

y(:,i) = y(:,i-1) + [(alpham*(1+(ylag1/p0)^n)-((a1)+(a2-a1)./(1+((r)./p1).^n1))*m)*timestep;
                     (alphap*(1/(1+(r/p2)^n2))*m - mup*p)*timestep;
                     (alphar*(1/(1+(p/p3)^n3))-murr*r)*timestep;
                     N*(alpham*Gx - M*Sx)*timestep + (N^0.5)*(timestep^0.5)*real((Sx*M + alpham*Gx)^0.5)*rnds(1);
                     N*(Fx*alphap*M - mup*P)*timestep + (N^0.5)*(timestep^0.5)*real((mup*P + alphap*Fx*M)^0.5)*rnds(2);
                     N*(alphar*Grx - murr*R)*timestep + (N^0.5)*(timestep^0.5)*real((murr*R + alphar*Grx)^0.5)*rnds(3)];
                                                            
i = i+1;
end

t=timestep:timestep:Tend;
pTraj1 = y(5,:)/N;
ix1 = find(pTraj1>20,1,'last');
Dist1(j) = t(ix1);

pTraj2 = y(2,:);
ix2 = find(pTraj2>20,1,'last');
Dist2(j) = t(ix2);

end

end

