function [ydot, datOut] = f(t,y,stp)

p = stp.p;
checkFigure = stp.checkFigure;

if size(y,2) == 1
    y = y';
end
ydot = zeros(size(y));
% -------- SIGNALING MODEL -----------

% b-AR module
LR = y(:,1).*y(:,2)./p(4);
LRG = LR.*y(:,3)./p(5);
RG = y(:,2).*y(:,3)./p(6);
BARKDESENS = p(7).*(LR+LRG);
BARKRESENS = p(8).*y(:,5);
PKADESENS = p(9).*y(:,17).*y(:,4);  
PKARESENS = p(10).*y(:,6);
GACT = p(11).*(RG+LRG);
HYD = p(12).*y(:,7);
REASSOC = p(13).*y(:,8).*y(:,9);
ydot(:,1) = p(1)-LR-LRG-y(:,1);
ydot(:,2) = y(:,4)-LR-LRG-RG-y(:,2);
ydot(:,3) = p(3)-LRG-RG-y(:,3);
ydot(:,4) = (BARKRESENS-BARKDESENS)+(PKARESENS-PKADESENS);
ydot(:,5) = BARKDESENS-BARKRESENS;
ydot(:,6) = PKADESENS-PKARESENS;
ydot(:,7) = GACT-HYD;
ydot(:,8) = HYD-REASSOC;
ydot(:,9) = GACT-REASSOC;
% end b-AR module

% cAMP module
Gsa_gtp_AC = y(:,10).*y(:,12)./p(27);
Fsk_AC = y(:,11).*y(:,12)./p(28);
AC_ACT_BASAL = p(19).*y(:,12).*p(15)./(p(23)+p(15));	    
AC_ACT_GSA = p(20).*Gsa_gtp_AC.*p(15)./(p(24)+p(15)); 
AC_ACT_FSK = p(21).*Fsk_AC.*p(15)./(p(25)+p(15));	   
PDE_ACT = p(22).*y(:,13).*y(:,16)./(p(26)+y(:,16));	
PDE_IBMX = y(:,13).*y(:,14)./p(29);
ydot(:,10) = y(:,7)-Gsa_gtp_AC-y(:,10);
ydot(:,11) = p(18)-Fsk_AC-y(:,11);
ydot(:,12) = p(14)-Gsa_gtp_AC-y(:,12);  
ydot(:,13) = p(16)-PDE_IBMX-y(:,13);
ydot(:,14) = p(17)-PDE_IBMX-y(:,14);
ydot(:,15) = AC_ACT_BASAL+AC_ACT_GSA+AC_ACT_FSK-PDE_ACT;
% end cAMP module

% PKA module
% 4/13/04 note: 
% I am currently not using an analytical solution for cAMPfree, but the
% below analytical solution does work.
% cAMPtemp = y(:,15)-(2.*A2RC_I+2.*A2R_I)-(2.*A2RC_II+2.*A2R_II);
% ydot(:,16) = cAMPtemp-sqrt(cAMPtemp.^2-4.*p(33).*(A2RC_I+A2RC_II))-y(:,16);
PKI = p(32).*p(36)./(p(36)+y(:,17)+y(:,18));
A2RC_I = (y(:,17)./p(35)).*y(:,17).*(1+PKI./p(36));
A2R_I = y(:,17).*(1+PKI./p(36));
A2RC_II = (y(:,18)./p(35)).*y(:,18).*(1+PKI./p(36));
A2R_II = y(:,18).*(1+PKI./p(36));
ARC_I = (p(33)./y(:,16)).*A2RC_I;
ARC_II = (p(33)./y(:,16)).*A2RC_II;
ydot(:,16) = y(:,15)-(ARC_I+2.*A2RC_I+2.*A2R_I)-(ARC_II+2.*A2RC_II+2.*A2R_II)-y(:,16);
PKAtemp = p(33).*p(34)./p(35)+p(33).*y(:,16)./p(35)+y(:,16).^2./p(35);
ydot(:,17) = 2.*p(30).*y(:,16).^2-y(:,17).*(1+PKI./p(36)).*(PKAtemp.*y(:,17)+y(:,16).^2);
ydot(:,18) = 2.*p(31).*y(:,16).^2-y(:,18).*(1+PKI./p(36)).*(PKAtemp.*y(:,18)+y(:,16).^2);
% end PKA module

% override d(cAMP_tot): hold it constant for Fig 2c
if strcmp(checkFigure, '2cd') || (strcmp(checkFigure, '2c') || strcmp(checkFigure, '2d') )
    ydot(:,15) = 0;
end


% PLB module
PLB = p(38)-y(:,19);
PLB_PHOSPH = p(41).*y(:,17).*PLB./(p(42)+PLB);
PLB_DEPHOSPH = p(43).*y(:,22).*y(:,19)./(p(44)+y(:,19));
ydot(:,19) = PLB_PHOSPH-PLB_DEPHOSPH;
 
Inhib1 = p(40)-y(:,20);
Inhib1p_PP1 = y(:,21).*y(:,22)./p(49);
Inhib1_PHOSPH = p(45).*y(:,17).*Inhib1./(p(46)+Inhib1); 
Inhib1_DEPHOSPH = p(47).*y(:,20)./(p(48)+y(:,20));
ydot(:,20) = Inhib1_PHOSPH-Inhib1_DEPHOSPH;
ydot(:,21) = y(:,20)-Inhib1p_PP1-y(:,21);
ydot(:,22) = p(39)-Inhib1p_PP1-y(:,22);

fracPLBp = y(:,19)./p(38);
fracPLB = PLB./p(38);
fracPLBo = stp.fracPLBo; % 0.9613; % adjust when changes are made to signaling model!
if strcmp(stp.fig6, 'LCC')
    ydot(:,19) = 0;
end
if stp.PLBKO
    fracPLB = zeros(size(fracPLB));
    fracPLBo = ones(size(fracPLBo));
end
% end PLB module

% LCC module
LCCa = p(50)-y(:,23);
LCCa_PHOSPH = p(37).*p(53).*y(:,18).*LCCa./(p(54) + p(37).*LCCa);
LCCa_DEPHOSPH = p(37).*p(57).*p(52).*y(:,23)./(p(58)+p(37).*y(:,23));
ydot(:,23) = LCCa_PHOSPH - LCCa_DEPHOSPH;
fracLCCap = y(:,23)./p(50);
fracLCCapo = 0.2041;
 
LCCb = p(50)-y(:,24);
LCCb_PHOSPH = p(37).*p(53).*y(:,18).*LCCb./(p(54)+p(37).*LCCb);   
LCCb_DEPHOSPH = p(37).*p(55).*p(51).*y(:,24)./(p(56)+p(37).*y(:,24));
ydot(:,24) = LCCb_PHOSPH-LCCb_DEPHOSPH;
fracLCCbp = y(:,24)./p(50);
fracLCCbpo = 0.2336;
if strcmp(stp.fig6,'PLB')
    ydot(:,23) = 0;
    ydot(:,24) = 0;
end
% end LCC module
% -------- END SIGNALING MODEL ---------
% -------- EC COUPLING MODEL -----------

% Constants
R = 8314;   % R     [J./kmol.*K]
Frdy = 96485;  % Frdy     [C./mol]
FoRT = Frdy./R./p(63);
zna = 1;    % Na valence
zk = 1;     % K valence
zca = 2;    % Ca valence
% Nernst Potentials
ena = (1./FoRT./zna).*log(p(64)./y(:,40));       % should be 70.54 mV
ek = (1./FoRT./zk).*log(p(65)./y(:,41));		 % should be -87.94 mV
eca = (1./FoRT./zca).*log(p(66)./y(:,42));  % should be 120 mV
%eks = (1./FoRT).*logn((ko+prnak.*nao)./(ki+prnak.*nai))	% should be -77.54
ecl = -40.0;

% I_Na: Fast Na Current
am = 0.32.*(y(:,43)+47.13)./(1-exp(-0.1.*(y(:,43)+47.13)));
bm = 0.08.*exp(-y(:,43)./11);
if y(:,43) >= -40
    ah = 0; aj = 0;
    bh = 1./(0.13.*(1+exp(-(y(:,43)+10.66)./11.1)));
    bj = 0.3.*exp(-2.535e-7.*y(:,43))./(1+exp(-0.1.*(y(:,43)+32)));
else
    ah = 0.135.*exp((80+y(:,43))./-6.8);
    bh = 3.56.*exp(0.079.*y(:,43))+3.1e5.*exp(0.35.*y(:,43));
    aj = (-1.2714e5.*exp(0.2444.*y(:,43))-3.474e-5.*exp(-0.04391.*y(:,43))).*(y(:,43)+37.78)./(1+exp(0.311.*(y(:,43)+79.23)));
    bj = 0.1212.*exp(-0.01052.*y(:,43))./(1+exp(-0.1378.*(y(:,43)+40.14)));
end
ydot(:,25) = 1e3.*(am.*(1-y(:,25))-bm.*y(:,25));
ydot(:,26) = 1e3.*(ah.*(1-y(:,26))-bh.*y(:,26));
ydot(:,27) = 1e3.*(aj.*(1-y(:,27))-bj.*y(:,27));
I_Na = p(67).*y(:,25).^3.*y(:,26).*y(:,27).*(y(:,43)-ena);

% I_Ca: L-type Calcium Current
alcc = 400.*exp( (y(:,43)+2)./10 );    % [1/sec]
blcc = 50.*exp( -1.*(y(:,43)+2)./13 );  % [1/sec]
flcc = p(72).*(0.375.*fracLCCap./fracLCCapo+0.625); % PHOSPHOREGULATION
ylccinf = 1.0./(1+exp((y(:,43)+55.0)./7.5 )) +0.1./(1+exp((-y(:,43)+21.0)./6.0 ));
tauylcc = 0.02 + 0.3./( 1 +exp( (y(:,43)+30)./9.5 ) );    % [sec]			
gamma = p(74).*y(:,42); % [1/sec./mM]		
vgamma = gamma.*((1-y(:,28)).^4+2.*y(:,28).*(1-y(:,28)).^3+4.*y(:,28).^2.*(1-y(:,28)).^2+8.*y(:,28).^3.*(1-y(:,28))+16.*y(:,28).^4.*(1-flcc./p(73)));
vomega = p(75).*((1-y(:,29)).^4+1./2.*y(:,29).*(1-y(:,29)).^3+1./4.*y(:,29).^2.*(1-y(:,29)).^2+1./8.*y(:,29).^3.*(1-y(:,29))+1./16.*y(:,29).^4);

ydot(:,28) = alcc.*(1-y(:,28))-blcc.*y(:,28);
ydot(:,29) = 2.*alcc.*(1-y(:,29))-blcc./2.*y(:,29);
ydot(:,30) = flcc.*(1-y(:,30))-p(73).*y(:,30);
ydot(:,31) = (ylccinf - y(:,31))./tauylcc;
ydot(:,32) = vomega.*(1-y(:,32))-vgamma.*y(:,32);

if abs(y(:,43)) < 1e-6
    disp('Warning! Voltage near zero, could influence ibarca and ibark!');
end
ibarca = p(76).*4.*(y(:,43).*Frdy.*FoRT) .* (1e-3.*exp(2.*y(:,43).*FoRT)-0.341.*p(66)) ./(exp(2.*y(:,43).*FoRT)-1);
ibark = p(77).*(y(:,43).*Frdy.*FoRT).*(y(:,41).*exp(y(:,43).*FoRT)-p(65)) ./(exp(y(:,43).*FoRT)-1);
favail = 0.5.*(0.4.*fracLCCbp./fracLCCbpo+0.60);   % PHOSPHOREGULATION
I_Ca = ibarca.*p(78).*favail.*y(:,28).^4.*y(:,30).*y(:,31).*y(:,32);
I_CaK = ibark./(1+I_Ca./p(79)).*p(78).*favail.*y(:,28).^4.*y(:,30).*y(:,31).*y(:,32);
%I_Ca = 0; I_CaK = 0;
I_Catot = I_Ca+I_CaK;


% I_to: Transient Outward K Current
rtoss = 1./(1+exp((y(:,43)+10.6)./-11.42));
stoss = 1./(1+exp((y(:,43)+45.3)./6.8841));
taurto = 1.0./(45.16.*exp(0.03577.*(y(:,43)+50))+98.9.*exp(-0.1.*(y(:,43)+38))); % [sec]
tausto = 0.35.*exp(-((y(:,43)+70)./15).^2)+0.035;	    % [sec]		 
taussto = 3.7.*exp(-((y(:,43)+70)./30).^2)+0.035;        % [sec]			 
ydot(:,33) = (rtoss-y(:,33))./taurto;
ydot(:,34) = (stoss-y(:,34))./tausto;
ydot(:,35) = (stoss-y(:,35))./taussto;
I_to = p(68).*y(:,33).*(0.886.*y(:,34)+0.114.*y(:,35)).*(y(:,43)-ek);   % [uA./uF]

% I_ss: Steady-state K Current
rssinf = 1./(1+exp(-(y(:,43)+11.5)./11.82));
taurss = 10./(45.16.*exp(0.03577.*(y(:,43)+50))+98.9.*exp(-0.1.*(y(:,43)+38)));
sssinf = 1./(1+exp((y(:,43)+87.5)./10.3));
tausss = 2.1;
ydot(:,36) = (rssinf-y(:,36))./taurss;
ydot(:,37) = (sssinf-y(:,37))./tausss;
I_ss = p(69).*y(:,36).*y(:,37).*(y(:,43)-ek);

% I_ki: Time-Independent K Current
aki = 1.02./(1+exp(0.2385.*(y(:,43)-ek-59.215)));
bki =(0.49124.*exp(0.08032.*(y(:,43)+5.476-ek)) + exp(0.06175.*(y(:,43)-ek-594.31))) ./(1 + exp(-0.5143.*(y(:,43)-ek+4.753)));
kiss = aki./(aki+bki);
I_ki = p(70).*sqrt(p(65)./5.4).*kiss.*(y(:,43)-ek) ;

% I_kp: Plateau K Current
kp = 1./(1+exp((7.488-y(:,43))./5.98));
I_kp = p(71).*kp.*(y(:,43)-ek);

% I_ncx: Na./Ca Exchanger Current
s4 = exp(p(84).*y(:,43).*FoRT).*y(:,40).^3.*p(66);
s5 = exp((p(84)-1).*y(:,43).*FoRT).*p(64).^3.*y(:,42);
I_ncx = p(80)./(p(81).^3+p(64).^3) ./(p(82)+p(66)) ./(1+p(83).*exp((p(84)-1).*y(:,43).*FoRT)) .*(s4-s5);

% I_nak: Na./K Pump Current
sigma = (exp(p(64)./67.3)-1)./7;
fnak = 1./(1+0.1245.*exp(-0.1.*y(:,43).*FoRT)+0.0365.*sigma.*exp(-y(:,43).*FoRT));
I_nak = p(85) .*fnak ./(1+(p(86)./y(:,40)).^1.5) .*(p(65)./(p(65)+p(87)));

% I_pca: Sarcolemmal Ca Pump Current
I_pca = p(88).*y(:,42)./(p(89)+y(:,42));
 
% I_cab: Ca Background Current
I_cab = p(90).*(y(:,43)-eca);
 
% I_nab: Na Background Current
I_nab = p(91).*(y(:,43)-ena);

% Total Membrane Currents
I_Na_tot = I_Na+I_nab+3.*I_ncx+3.*I_nak;          % [uA/uF]
I_K_tot = I_to+I_ss+I_ki+I_kp-2.*I_nak+I_CaK;    % [uA/uF]
I_Ca_tot = I_Ca+I_cab+I_pca-2.*I_ncx;            % [uA/uF]

% Calcium Induced Calcium Release (CICR)
trel = y(:,44)+2e-3;
ryron = 1-exp(-trel./p(97));  
ryroff = exp(-trel./p(98));
grel = p(99)./(1+exp((I_Ca_tot+5)./0.9));         % adjusted for rat  [1/sec] 
I_rel = grel.*ryron.*ryroff.*(y(:,39)-y(:,42));        %   [mM/sec]

% Other SR fluxes and concentrations
Km_up = p(95).*(1+2.*fracPLB)./(1+2.*fracPLBo);     % PHOSPHOREGULATION
I_up = p(94).*y(:,42).^2./(Km_up.^2+y(:,42).^2);         %   [mM/sec]
% Original: I_up = p(94).*y(:,42).^2./(p(95).^2+y(:,42).^2);         %   [mM/sec]
I_leak = p(94).*y(:,38)./p(96);                     %   [mM/sec]
I_tr = (y(:,38)-y(:,39))./p(105);                    %   [mM/sec]
Bjsr = 1./( 1+p(103).*p(104)./(p(104)+y(:,39)).^2 );
ydot(:,38) = I_up-I_leak-I_tr.*p(61)./p(60);        %   [mM/sec]
ydot(:,39) = Bjsr.*(I_tr-I_rel);                   %   [mM/sec]
SRcontent = 1e3.*((y(:,39)+y(:,39)./Bjsr).*p(61)./p(59)+y(:,38).*p(60)./p(59));    % [umol/L cytosol]

% Cytoplasmic Calcium Buffering
btrpn = p(106).*p(109)./(p(109)+y(:,42)).^2;
bcmdn = p(107).*p(110)./(p(110)+y(:,42)).^2;
bindo = p(108).*p(111)./(p(111)+y(:,42)).^2;
Bmyo = 1./( 1+ bcmdn + btrpn + btrpn + bindo);

% Ion Concentrations and Membrane Potential
ydot(:,40) = -1e3.*I_Na_tot.*p(62)./(p(59).*zna.*Frdy);          % [mM/sec] 
ydot(:,41) = -1e3.*I_K_tot.*p(62)./(p(59).*zk.*Frdy);            % [mM/sec]
ydot(:,42) = -Bmyo.*(1e3.*I_Ca_tot.*p(62)./(p(59).*zca.*Frdy) ...
    +((I_up-I_leak).*p(60)./p(59))-(I_rel.*p(61)./p(59)));    % [mM/sec]

% Simulation type
% protocol = 'pace';

protocol = stp.protocol;
I_app = zeros(size(t));
switch protocol
    case 0 % none
%         I_app = zeros(size(t));
    case 1 %'pace',        % pace w./ current injection at rate 'rate'
		rate = 1;
        for i=1:length(t)
            if mod(t(i)+0.9,1./rate) <= 5e-3
                I_app(i) = 10.0;
            else
                I_app(i) = 0.0;
            end
        end
    case 2 %'vclamp',      
		V_hold = -55;
        V_test = -10;
        for i=1:length(t)
            if (t(i) > 59.1 & t(i) < 59.5)
                V_clamp = V_test;
            else
                V_clamp = V_hold;
            end
            R_clamp = .02; 
            I_app(i) = (V_clamp-y(i,43))./R_clamp;
        end
end  

ydot(:,43) = -1e3.*(I_Ca_tot+I_K_tot+I_Na_tot-I_app);

% CICR timing: y(:,44) tracks the last time when Vdot > 30 mV./msec
if (ydot(:,43) > 30e3) 
    ydot(:,44) = 1-1e4.*y(:,44); 
else
    ydot(:,44) = 1;
end

% ----- END EC COUPLING MODEL ---------------


% Shelley: output currents in ydot
% I_Na+I_nab+3.*I_ncx+3.*I_nak;          % [uA./uF]
% I_to+I_ss+I_ki+I_kp-2.*I_nak+I_CaK;    % [uA./uF]
% I_Ca+I_cab+I_pca-2.*I_ncx;            % [uA./uF]
% I_up-I_leak

if min(size(ydot)) == 1
    ydot = ydot';
end

datOut = struct;
datOut.I_Ca = I_Ca;
datOut.I_app = I_app;
datOut.PKA_DESENS = PKADESENS; 
datOut.PKA_RESENS = PKARESENS;
datOut.PKI = PKI;
datOut.fracPLBp = fracPLBp;
datOut.fracPLB = fracPLB;
datOut.Gsa_gtp_AC = Gsa_gtp_AC;