%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = ds_dt(t,y,options,eNa,eK,gNa,gK)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

V = y(:,:,1);
Ca = y(:,:,2);
x = y(:,:,3);
m = y(:,:,4);
h = y(:,:,5);
j = y(:,:,6);
d = y(:,:,7);
f = y(:,:,8);

%gK% = 0.35;
%gNa% = 4.0;

%eK% = -85;
%eNa% = 50;
dy = 0*y;
%gates
c4 = 0; c5 = 0;
ax = (0.0005*exp(0.0830*(V+50)) + c4*(V+c5))./(exp(0.057*(V+50))+1);
bx = (0.0013*exp(-0.060*(V+20)) + c4*(V+c5))./(exp(-0.04*(V+20))+1);
dy(:,:,3) = ax.*(1-x)-bx.*x;
am = (0.0000*exp(0.0000*(V+47)) -  1*(V+47))./(exp(-0.1*(V+47))-1);
bm = (40.000*exp(-0.056*(V+72)) + c4*(V+c5))./(exp(0.00*(V+72))+0);
dy(:,:,4) = am.*(1-m)-bm.*m;
ah = (0.1260*exp(-0.250*(V+77)) + c4*(V+c5))./(exp(0.0*(V+77))+0);
bh = (1.7000*exp(0.0000*(V+22.5))+c4*(V+c5))./(exp(-0.082*(V+22.5))+1);
dy(:,:,5) = ah.*(1-h)-bh.*h;
aj = (0.0550*exp(-0.250*(V+78)) + c4*(V+c5))./(exp(-0.2*(V+78))+1);
bj = (0.3000*exp(0.0000*(V+32)) + c4*(V+c5))./(exp(-0.1*(V+32))+1);
dy(:,:,6) = aj.*(1-j)-bj.*j;
ad = (0.0950*exp(-0.010*(V-5))  + c4*(V+c5))./(exp(-0.072*(V-5))+1);
bd = (0.0700*exp(-0.017*(V+44)) + c4*(V+c5))./(exp(0.05*(V+44))+1);
dy(:,:,7) = ad.*(1-d)-bd.*d;
af = (0.0120*exp(-0.008*(V+28)) + c4*(V+c5))./(exp(0.15*(V+28))+1);
bf = (0.0065*exp(-0.020*(V+30)) + c4*(V+c5))./(exp(-0.2*(V+30))+1);
dy(:,:,8) = af.*(1-f)-bf.*f;

%currents
iK = gK*(4*(exp(0.04*(V-eK))-1)./(exp(0.08*(V+53))+exp(0.04*(V+53)))+0.2*(V+23)./(1-exp(-0.04*(V+23))));
iX = x*0.8.*(exp(0.04*(V+77))-1)./exp(0.04*(V+35));
gNaC = 0.003;
iNa = (gNa*m.^3.*h.*j+gNaC).*(V-eNa);
eS = -82.3-13.0287*log(Ca); gS = 0.09;
iS = gS*d.*f.*(V-eS);

dy(:,:,2) = -1e-7*iS+0.07*(1e-7-Ca);
dy(:,:,1) = -(1/1)*(iK+iX+iNa+iS);

