- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-09-03 12:15:33+12:00
- Desc:
- Init
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/6fa/rawfile/e459ccc5beeea2965725c7eea6409c97c3bbd563/matlab/f.m
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;
