function [ y, dy] = ChemicalLangevinEquationFinal( tau,v,Tend,timestep,parms )
%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;

y = zeros(6,ttot);
dy = zeros(6,ttot);
hist = repmat(v,1,h);
y(:,[1:h])=hist;

for i =((h+1):ttot)

% dy=zeros(9,1);

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

G = 1/(1+(ylag1/p0)^n);
Gr = 1/(1+(p/p3)^n3);
F = 1/(1+(r/p2)^n2);

% a1=a1;a2=a2;
A=(a2)-(a1);
S=(a1)+A./(1+((r)./p1).^n1);

% mrna
dy(1,i) = alpham*G - m*S;

% protein
dy(2,i) = F*alphap*m - mup*p; 

% mir-9
dy(3,i) = alphar*Gr - murr*r;

% Gd =  -alpham*(n/ylag1)*((ylag1/p0)^n)/(1+(ylag1/p0)^n)^2;
% Sd = -(a2-a1)*(n1/r)*((r/p1)^n1)/(1+(r/p1)^n1)^2;
% Fd =  -alphap*(n2/r)*((r/p2)^n2)/(1+(r/p2)^n2)^2;
% Grd =  -alphar*(n3/p)*((p/p3)^n3)/(1+(p/p3)^n3)^2;
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);

% a1=a1;a2=a2;
A=(a2)-(a1);
Sx=(a1)+A./(1+((R)./p1).^n1);

dy(4,i) = N*(alpham*Gx - M*Sx);
dy(5,i) = N*(Fx*alphap*M - mup*P); 
dy(6,i) = N*(alphar*Grx - murr*R);

etam2 = real((S*M + alpham*Gx)^0.5);
etap2 = real((mup*P + alphap*Fx*M)^0.5);
etar2 = real((murr*R + alphar*Grx)^0.5);

etam = normrnd(0,etam2);
etap = normrnd(0,etap2);
etar = normrnd(0,etar2);

% dy(4,i) = N*dy(1,i) + (N^0.5)*dy(4,i);
% dy(5,i) = N*dy(2,i) + (N^0.5)*dy(5,i);
% dy(6,i) = N*dy(3,i) + (N^0.5)*dy(6,i);

y(1,i) = y(1,(i-1))+dy(1,i)*timestep;
y(2,i) = y(2,(i-1))+dy(2,i)*timestep;
y(3,i) = y(3,(i-1))+dy(3,i)*timestep;
y(4,i) = y(4,(i-1))+dy(4,i)*timestep + (N^0.5)*(timestep^0.5)*etam;
y(5,i) = y(5,(i-1))+dy(5,i)*timestep+ (N^0.5)*(timestep^0.5)*etap;
y(6,i) = y(6,(i-1))+dy(6,i)*timestep+ (N^0.5)*(timestep^0.5)*etar;

dy(4,i) = dy(4,i) + (N^0.5)*etam;
dy(5,i) = dy(4,i) + (N^0.5)*etam;
% dy(6,i) = dy(4,i) + (N^0.5)*etam;
dy(6,i) = N*(alphar*Grx);
dy(3,i) = alphar*Gr;
end
end

