INa and IKv4.3 heterogeneity in canine LV myocytes (Flaim et al 2006)

 Download zip file 
Help downloading and running models
Accession:83491
"The roles of sustained components of INa and IKv43 in shaping the action potentials (AP) of myocytes isolated from the canine left ventricle (LV) have not been studied in detail. Here we investigate the hypothesis that these two currents can contribute substantially to heterogeneity of early repolarization and arrhythmic risk.... The resulting simulations illustrate ways in which KChIP2- and Ca2+- dependent control of IKv43 can result in a sustained outward current that can neutralize INaL in a rate- and myocyte subtype-dependent manner. Both these currents appear to play significant roles in modulating AP duration and rate dependence in midmyocardial myocytes. ... By design, these models allow upward integration into organ models or may be used as a basis for further investigations into cellular heterogeneities." See paper for more and details.
Reference:
1 . Flaim SN, Giles WR, McCulloch AD (2006) Contributions of sustained INa and IKv43 to transmural heterogeneity of early repolarization and arrhythmogenesis in canine left ventricular myocytes. Am J Physiol Heart Circ Physiol 291:H2617-29 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Electrogenic pump;
Brain Region(s)/Organism:
Cell Type(s): Heart cell; Cardiac ventricular cell;
Channel(s): I Na,t; I A; I K; Late Na; Na/Ca exchanger; I_Na,Ca; I_SERCA; Na/K pump;
Gap Junctions:
Receptor(s):
Gene(s): Kv4.3 KCND3;
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Ion Channel Kinetics; Action Potentials; Heart disease; Sodium pump; Markov-type model;
Implementer(s): Flaim, Sarah [flaim at comlab.ox.ac.uk];
Search NeuronDB for information about:  I Na,t; I A; I K; Late Na; Na/Ca exchanger; I_Na,Ca; I_SERCA; Na/K pump;
/
flaim_ajp
NewerICs
readme.txt
fcn_fgm.m
Script_FGM.m
                            
%        -------------------------------------------------------------
%
%          NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE
%
%          Copyright 2005, The Johns Hopkins University
%             School of Medicine and Whiting School of Engineering.
%             All rights reserved.
%             For research use only; commercial use prohibited.
%             Distribution without permission of Joseph L. Greenstein or
%             Raimond L. Winslow not permitted.
%           (jgreenst@jhu.edu, rwinslow@jhu.edu)
%
%          Copyright 2006,  University of California, San Diego
%            Department of Bioengineering
%            All rights reserved.
%            For research use only; commercial use prohibited.
%            Distribution without permission of Sarah N. Flaim or
%            Andrew D. McCulloch not permitted.
%           (shealy@bioeng.ucsd.edu, amcculloch@ucsd.edu)      
%
%          Original model and manuscript:
%          Name of Program: 40-State coupled LCC-RyR Model
%          Version: version 1.0
%          Date: September 2005
%         
%          Greenstein, J. L., R. Hinch, and R.L. Winslow. 2006. Protein Geometry and Placement in the
%          Cardiac Dyad Influence Macroscopic Properties of Calcium-Induced Calcium-Release. 
%          Biophys J. 90(1): 77-91.
%
%          Modified model and manuscript:
%          Description: Canine ventricular myocyte models of endo, epi and
%          M cell ionic current and excitation-contraction coupling.
%          Major modifications: Replacement of Hodgkin-Huxley INa model with Clancy Markov model, including late INa
%                              Inclusion of novel open state inactivation for Ito1 model
%                              Inclusion of selected transmurally varying ion channel parameters                                       
%          Version: version 2.0
%          Date: March 2006
%        
%          Flaim, S.N, Giles, W.R., and McCulloch, A.D 2006. Contributions of sustained INa and IKv4.3 to 
%          transmural heterogeneity of early repolarization and arrhythmogenesis in canine left ventricular 
%          myocytes
%          Am J Physiol Heart Circ Physiol 291: H2617-H2629, 2006
%
%
%        ------------------------------------------------------------- 
  

function dy = fcn_fgm(t,y)

global Vclamp_flag Period Shift Pulse_duration Pulse_amplitude Nstates_CaRU IKs_array
global Counter T_array ICaL_array JCaL_array JRyR_array CaSSavg_array celltype

NCaRU = 50000;                  %   Number of Ca Release Units
Faraday = 96.5;                 %   Faraday's constant (C/mmol)
Temp = 310;                     %   Temperature (K) 310 = 37, 296 = 23
Rgas = 8.314;                   %   Universal Gas constant (J/[mol K])
RT_over_F = Rgas*Temp/Faraday;
Csa = 153.4e-6;                 %   Membrane capacitance (pF)
VSS = 0.203e-12;                %   Single dyad subspace volume (microL)
VSR = 1.3913e-6;                %   SR volume (microL)
Vmyo = 25.84e-6;                %   Myoplasmic volume (microL)
PCa = 9.13e-13;                 %   LCC single channel permeability (cm^3/s)
JRyRmax = 5 * 3.92;             %   Ca flux rate through open RyR (1/ms)
r_xfer = 220;                   %   Ca flux rate from SS to myoplasm (1/ms)
Cao = 2;                        %   Extracellular [Ca] (mM)
Ko = 4;                         %   Extracellular [K] (mM)
Nao = 138;                      %   Extracellular [Na] (mM)
Clo = 150;                      %   Extracellular [Cl] (mM)
Cli = 20;                       %   Intracellular [Cl] (mM)

KdIto2 = 0.1502;                %   Dissociation constant of Cl channel (Ito2) for Ca (mM)
PCl = 2.65e-15;                 %   Cl single channel permeability (cm^3/s)

i768v = 1.0;                    % Mutation = 1.0 for normal, 1.5 for mutation.    
switch celltype
    case 'mid'                       
        GKs = 0.01;                    %    Maximal conductance of IKs (mS/microF)    
        kNaCa = 0.27;                  %    Scale factor for Na/Ca exchanger (A/F) - potentially heterogeneous
        KSR = 1.0;                     %    Scale factor for SR Ca ATPase
        GNal = 2.3;%2.3;               %    Maximal conductance of INal(mS/uF). New parameter.
        Kv4p3Scale = 2.5;              %    Scale factor for kv4.3 current. New parameter.
        betaa0Kv43 = 0.080185;         %    Scale factor for kv4.3 current. New parameter.
        Ito2scale = 0.5;               %    Scale factor for icacl current. New parameter.
        k_alphai = 1.74;               %    Scale factor for kv4.3 current. New parameter.
        k_betai = 1.41;                %    Scale factor for kv4.3 current. New parameter.
        k_alphaoi = 750;               %    Scale factor for kv4.3 current. New parameter.
        k_betaoi = 12*5.2310e-4;       %    Scale factor for kv4.3 current. New parameter.
        betaoi = 0.038;                %    Scale factor for kv4.3 current. New parameter.
        Cm = 153;                      %    Capacitance - not used. 
    case 'epi'
        GKs = 0.02;                    %    Maximal conductance of IKs (mS/microF)   
        kNaCa = 0.27;                  %    Scale factor for Na/Ca exchanger (A/F) - potentially heterogeneous
        KSR = 2;%1.0*2.0;              %    Scale factor for SR Ca ATPase         
        GNal = 1.15;                   %    Maximal conductance of INal(mS/uF) 
        Kv4p3Scale = 3.0;              %    Scale factor for kv4.3 current. New parameter.
        k_alphai = 2.25;               %    Scale factor for kv4.3 current. New parameter.         
        k_betai = 1.125;               %    Scale factor for kv4.3 current. New parameter.
        k_alphaoi = 1000;              %    Scale factor for kv4.3 current. New parameter.
        k_betaoi = 16*5.2310e-4;       %    Scale factor for kv4.3 current. New parameter.
        betaoi = 0.038;                %    Scale factor for kv4.3 current. New parameter.
        betaa0Kv43 = 0.080185;         %    Scale factor for kv4.3 current. New parameter.
        Ito2scale = 0.5;               %    Scale factor for icacl current. New parameter.
        Cm = 182;                      %   Capacitance - not used. 
    case 'endo'
        GKs = 0.02;                    %   Maximal conductance of IKs (mS/microF)    
        kNaCa = 0.27;                  %   Scale factor for Na/Ca exchanger (A/F)
        KSR = 1.0;                     %   Scale factor for SR Ca ATPase
        GNal = 1.15;                   % Maximal conductance of INal(mS/uF) 
        Kv4p3Scale = 0.5;              %   Scale factor for kv4.3 current. New parameter.
        k_alphai = 0.18;               %   Scale factor for kv4.3 current. New parameter.
        k_betai = 2.25;                %   Scale factor for kv4.3 current. New parameter.
        k_alphaoi = 500;               %   Scale factor for kv4.3 current. New parameter.
        k_betaoi = 8*5.2310e-4;        %   Scale factor for kv4.3 current. New parameter.
        betaoi = 0.038;                %   Scale factor for kv4.3 current. New parameter.
        betaa0Kv43 = 0.080185;         %   Scale factor for kv4.3 current. New parameter.
        Ito2scale = 2.5;               %    Scale factor for icacl current. New parameter.
        Cm = 124;                      %   Capacitance - not used.
end

GKr = 0.029;                    %   Maximal conductance of IKr (mS/microF)
%GKs = 0.035;                    %   Maximal conductance of IKs (mS/microF)
GK1 = 2.73;                     %   Maximal conductance of IK1 (mS/microF)
GKp = 0.003;                    %   Maximal conductance of IKp (mS/microF)
GNa = 12.8;                     %   Maximal conductance of INa (mS/microF)
%kNaCa = 0.27;                   %   Scale factor for Na/Ca exchanger (A/F)
KmNa = 87.5;                    %   Na half saturation constant for Na/Ca exchanger (mM)
KmCa = 1.38;                    %   Ca half saturation constant for Na/Ca exchanger (mM)
KmK1 = 13.0;                    %   K half saturation constant for IK1 (mM)
ksat = 0.2;                     %   Na/Ca exchange saturation factor at negative potentials
eta = 0.35;                     %   Controls voltage dependence of Na/Ca exchange
INaKmax = 1.7325;               %   Scale factor for Na/K pump (A/F)
KmNai = 20;                     %   Na half saturation constant for Na/K pump (mM)
KmKo = 1.5;                     %   K half saturation constant for Na/K pump (mM)
IpCamax = 0.03;                 %   Maximum sarcolemmal Ca pump current (A/F)
KmpCa = 0.0005;                 %   Ca half saturation constant for sarcolemmal Ca pump (mM)
GCab = 0;                       %   Maximal conductance of ICab (mS/microF)
GNab = 0.00001;                 %   Maximal conductance of INab (mS/microF)

Kfb = 0.26e-3;                  %   Forward half-saturation constant for SR Ca ATPase (mM)
Krb = 1.8;                      %   Reverse half-saturation constant for SR Ca ATPase (mM)
%KSR = 1;                        %   Scale factor for SR Ca ATPase
Nfb = 0.75;                     %   Forward cooperativity constant for SR Ca ATPase
Nrb = 0.75;                     %   Reverse cooperativity constant for SR Ca ATPase
vmaxf = KSR*2.0961e-4;              %   SR Ca ATPase forward rate parameter (mM/ms)
vmaxr = KSR*2.0961e-4;              %   SR Ca ATPase reverse rate parameter (mM/ms)

KvScale = 1.5876;               %   Scale factor for Kv4.3 and Kv1.4 currents
Kv43Frac = 0.85;                %   Fraction of Ito1 (peak current at 40 mV) that is made up of Kv4.3 
%Kv43Frac = 0.70;                %   Fraction of Ito1 (peak current at 40 mV) that is made up of Kv4.3 
GKv43 = Kv43Frac*KvScale*0.1*Kv4p3Scale;   %   Maximal conductance of IKv43 (mS/microF)
PKv14 = (1-Kv43Frac)*KvScale*4.792933e-7;       %   Maximal conductance of IKv14 (mS/microF)


alphaa0Kv43 = 0.543708 ;         %   IKv43 rate parameter (1/ms)
aaKv43 = 0.028983;              %   IKv43 rate parameter (1/mV)
%betaa0Kv43 = 0.080185;          %   IKv43 rate parameter (1/ms)
baKv43 = 0.0468437/2;             %   IKv43 rate parameter (1/mV)
alphai0Kv43 = (0.0498424)*k_alphai;        %   IKv43 rate parameter (1/ms)
aiKv43 = 0.000373016;           %   IKv43 rate parameter (1/mV)
betai0Kv43 = (0.000819482)*k_betai;       %   IKv43 rate parameter (1/ms)
biKv43 = 0.00000005374;         %   IKv43 rate parameter (1/mV)
alphaa0Kv14 = 6*0.31551;        %   IKv14 rate parameter (1/ms)
aaKv14 = 0.00695;               %   IKv14 rate parameter (1/mV)
betaa0Kv14 = 6*0.001966;        %   IKv14 rate parameter (1/ms)
baKv14 = 0.08527;               %   IKv14 rate parameter (1/mV)
alphai0Kv14 = 6*0.0004938;      %   IKv14 rate parameter (1/ms)
betai0Kv14 = 6*0.0000176;       %   IKv14 rate parameter (1/ms)
f1Kv43 = 1.8936;                %   IKv43 rate parameter
f2Kv43 = 14.224647456;          %   IKv43 rate parameter
f3Kv43 = 158.574378389;         %   IKv43 rate parameter
f4Kv43 = 142.936645351;         %   IKv43 rate parameter
b1Kv43 = 6.77348;               %   IKv43 rate parameter
b2Kv43 = 15.6212705152;         %   IKv43 rate parameter
b3Kv43 = 28.7532603313;         %   IKv43 rate parameter
b4Kv43 = 524.576206679;         %   IKv43 rate parameter
f1Kv14 = 0.20005;               %   IKv14 rate parameter
f2Kv14 = 0.320280;              %   IKv14 rate parameter
f3Kv14 = 13.50909223;           %   IKv14 rate parameter
f4Kv14 = 1151.7651385;          %   IKv14 rate parameter
b1Kv14 = 2.2300;                %   IKv14 rate parameter
b2Kv14 = 12.000299;             %   IKv14 rate parameter
b3Kv14 = 5.3701338025;          %   IKv14 rate parameter
b4Kv14 = 5.2396395511;          %   IKv14 rate parameter

A0_HERG = 0.017147641733086;    %   IKr rate parameter (1/ms)
B0_HERG = 0.03304608038835;     %   IKr rate parameter (1/mV)
A1_HERG = 0.03969328381141;     %   IKr rate parameter (1/ms)
B1_HERG = -0.04306054163980;    %   IKr rate parameter (1/mV)
A2_HERG = 0.02057448605977;     %   IKr rate parameter (1/ms)
B2_HERG = 0.02617412715118;     %   IKr rate parameter (1/mV)
A3_HERG = 0.00134366604423;     %   IKr rate parameter (1/ms)
B3_HERG = -0.02691385498399;    %   IKr rate parameter (1/mV)
A4_HERG = 0.10666316491288;     %   IKr rate parameter (1/ms)
B4_HERG = 0.00568908859717;     %   IKr rate parameter (1/mV)
A5_HERG = 0.00646393910049;     %   IKr rate parameter (1/ms)
B5_HERG = -0.04536642959543;    %   IKr rate parameter (1/mV)
A6_HERG = 0.00008039374403;     %   IKr rate parameter (1/ms)
B6_HERG = 0.00000069808924;     %   IKr rate parameter (1/mV)
T_Const_HERG = 5.320000001;     %   IKr rate scale factor for temperature
C2H_to_C3H = T_Const_HERG*0.02608362043337;     %   IKr transition rate (1/ms)
C3H_to_C2H = T_Const_HERG*0.14832978132145;     %   IKr transition rate (1/ms)


CSQNtot = 2.7;                  %   Total SR calsequestrin concentration (mM)
KmCSQN = 0.63;                  %   Ca half saturation constant for calsequestrin (mM)
LTRPNtot = 0.07;                %   Total troponin low-affinity site concentration (mM)
HTRPNtot = 0.14;                %   Total troponin high-affinity site concentration (mM)
khtrpn_plus = 20;               %   Troponin Ca binding rate parameter (1/[mM ms])
khtrpn_minus = 0.066e-3;        %   Troponin Ca binding rate parameter (1/ms)
kltrpn_plus = 40;               %   Troponin Ca binding rate parameter (1/[mM ms])
kltrpn_minus = 0.04;            %   Troponin Ca binding rate parameter (1/ms)
CMDNtot = 0.05;                 %   Total myoplasmic calmodulin concentration (mM)
EGTAtot = 0;                    %   Total myoplasmic EGTA concentration (mM)
KmCMDN = 2.38e-3;               %   Ca half saturation constant for calmodulin (mM)
KmEGTA = 1.5e-4;                %   Ca half saturation constant for EGTA (mM)

V = y(Nstates_CaRU+1);          %   Assign states to local variables
P_UO = y(Nstates_CaRU+2); %mNa = y(Nstates_CaRU+2);
P_UIF = y(Nstates_CaRU+3); %hNa = y(Nstates_CaRU+3);
P_UC2 = y(Nstates_CaRU+4); %jNa = y(Nstates_CaRU+4);
Nai = y(Nstates_CaRU+5);
Ki = y(Nstates_CaRU+6);
Cai = y(Nstates_CaRU+7);
CaSR = y(Nstates_CaRU+8);
LTRPNCa = y(Nstates_CaRU+9);
HTRPNCa = y(Nstates_CaRU+10);
xKs = y(Nstates_CaRU+11);
C0Kv43 = y(Nstates_CaRU+12);
C1Kv43 = y(Nstates_CaRU+13);
C2Kv43 = y(Nstates_CaRU+14);
C3Kv43 = y(Nstates_CaRU+15);
OKv43 = y(Nstates_CaRU+16);
CI0Kv43 = y(Nstates_CaRU+17);
CI1Kv43 = y(Nstates_CaRU+18);
CI2Kv43 = y(Nstates_CaRU+19);
CI3Kv43 = y(Nstates_CaRU+20);
OIKv43 = y(Nstates_CaRU+21);
C0Kv14 = y(Nstates_CaRU+22);
C1Kv14 = y(Nstates_CaRU+23);
C2Kv14 = y(Nstates_CaRU+24);
C3Kv14 = y(Nstates_CaRU+25);
OKv14 = y(Nstates_CaRU+26);
CI0Kv14 = y(Nstates_CaRU+27);
CI1Kv14 = y(Nstates_CaRU+28);
CI2Kv14 = y(Nstates_CaRU+29);
CI3Kv14 = y(Nstates_CaRU+30);
OIKv14 = y(Nstates_CaRU+31);
C1Herg = y(Nstates_CaRU+32);
C2Herg = y(Nstates_CaRU+33);
C3Herg = y(Nstates_CaRU+34);
OHerg = y(Nstates_CaRU+35);
IHerg = y(Nstates_CaRU+36);
P_UC3 = y(Nstates_CaRU+37); %hNal = y(Nstates_CaRU+37);
OIKv43Ca = y(Nstates_CaRU+38);
P_UIM1 = y(Nstates_CaRU+39);
P_UC1 = y(Nstates_CaRU+40);
P_IC3 = y(Nstates_CaRU+41);
P_IC2 = y(Nstates_CaRU+42);
P_UIM2 = y(Nstates_CaRU+43);
P_LC3 = y(Nstates_CaRU+44);
P_LC2 = y(Nstates_CaRU+45);
P_LC1 = y(Nstates_CaRU+46);
P_LO = y(Nstates_CaRU+47);


VF_over_RT=V/RT_over_F;
VFsq_over_RT=(1000*Faraday)*VF_over_RT;

time_on_Is1 = floor((t-Shift)/Period)*Period;       %   Establish stimulus current
time_off_Is1 = time_on_Is1+Pulse_duration;
Istim = 0;
if ((t-Shift) >= time_on_Is1 & (t-Shift) <= time_off_Is1)
    Istim = Istim + Pulse_amplitude;
end

fL = 0.85;                      %   L-type Ca channel rate parameter (1/ms)
gL = 2;                         %   L-type Ca channel rate parameter (1/ms)
bL = 32.1948;                   %   L-type Ca channel rate parameter
aL = 12.8878;                   %   L-type Ca channel rate parameter
omega = 0.0269659;              %   L-type Ca channel rate parameter (1/ms)
gammacf = 1.39;                 %   L-type Ca channel rate parameter (1/[mM ms])

alpha = 0.835399*exp(0.0269241*(V-35));     %   L-type Ca channel rate parameter (1/ms)
beta =  0.0331584*exp(-0.0934594*(V-35));   %   L-type Ca channel rate parameter (1/ms)

yCainf = 0.4/(1+exp((V+12.5)/5)) + 0.6;     %   L-type Ca channel steady-state inactivation
tauyCa = 340/(1+exp((V+30)/12)) + 60;       %   L-type Ca channel inactivation time constant (ms)
kby = yCainf/tauyCa;            %   L-type Ca channel rate parameter (1/ms)
kfy = (1 - yCainf)/tauyCa;      %   L-type Ca channel rate parameter (1/ms)

k12 = 5 * 877.5;                %   RyR channel rate parameter (1/ms)
k21 = 5 * 250;                  %   RyR channel rate parameter (1/ms)
k23 = 2.358e8;                  %   RyR channel rate parameter (1/ms)
k32 = 9.6;                      %   RyR channel rate parameter (1/ms)
k34 = 1.415e6;                  %   RyR channel rate parameter (1/ms)
k43 = 13.65;                    %   RyR channel rate parameter (1/ms)
k45 = 0.07;                     %   RyR channel rate parameter (1/ms)
k54 = 93.385;                   %   RyR channel rate parameter (1/ms)
k56 = 1.887e7;                  %   RyR channel rate parameter (1/ms)
k65 = 30;                       %   RyR channel rate parameter (1/ms)
k25 = 2.358e6;                  %   RyR channel rate parameter (1/ms)
k52 = 0.001235;                 %   RyR channel rate parameter (1/ms)

CaSS_1 = Cai;                   %   Ca concentration for each open-closed configuaration
CaSS_2 = (PCa*2*VF_over_RT/VSS*0.341*Cao/(exp(2*VF_over_RT)-1) + r_xfer*Cai) ...
          / (PCa*2*VF_over_RT/VSS/(1-exp(-2*VF_over_RT)) + r_xfer);
CaSS_3 = (JRyRmax*CaSR + r_xfer*Cai)/(JRyRmax + r_xfer);
CaSS_4 = (PCa*2*VF_over_RT/VSS*0.341*Cao/(exp(2*VF_over_RT)-1) + JRyRmax*CaSR + r_xfer*Cai) ...
          / (PCa*2*VF_over_RT/VSS/(1-exp(-2*VF_over_RT)) + JRyRmax + r_xfer);

LCC_1_to_2 = alpha;             %   Set L-type Ca channel rates (1/ms)
LCC_2_to_1 = beta;
LCC_2_to_3 = fL;
LCC_3_to_2 = gL;
LCC_1_to_4_CaSS_1 = gammacf*CaSS_1;
LCC_1_to_4_CaSS_2 = gammacf*CaSS_2;
LCC_1_to_4_CaSS_3 = gammacf*CaSS_3;
LCC_1_to_4_CaSS_4 = gammacf*CaSS_4;
LCC_4_to_1 = omega;
LCC_4_to_5 = aL*alpha;
LCC_5_to_4 = beta/bL;
LCC_2_to_5_CaSS_1 = aL*gammacf*CaSS_1;
LCC_2_to_5_CaSS_2 = aL*gammacf*CaSS_2;
LCC_2_to_5_CaSS_3 = aL*gammacf*CaSS_3;
LCC_2_to_5_CaSS_4 = aL*gammacf*CaSS_4;
LCC_5_to_2 = omega/bL;

LCC_6_to_7 = alpha;
LCC_7_to_6 = beta;
LCC_7_to_8 = fL;
LCC_8_to_7 = gL;
LCC_6_to_9_CaSS_1 = gammacf*CaSS_1;
LCC_6_to_9_CaSS_2 = gammacf*CaSS_2;
LCC_6_to_9_CaSS_3 = gammacf*CaSS_3;
LCC_6_to_9_CaSS_4 = gammacf*CaSS_4;
LCC_9_to_6 = omega;
LCC_9_to_10 = aL*alpha;
LCC_10_to_9 = beta/bL;
LCC_7_to_10_CaSS_1 = aL*gammacf*CaSS_1;
LCC_7_to_10_CaSS_2 = aL*gammacf*CaSS_2;
LCC_7_to_10_CaSS_3 = aL*gammacf*CaSS_3;
LCC_7_to_10_CaSS_4 = aL*gammacf*CaSS_4;
LCC_10_to_7 = omega/bL;

LCC_1_to_6 = kfy;
LCC_2_to_7 = kfy;
LCC_3_to_8 = kfy;
LCC_4_to_9 = kfy;
LCC_5_to_10 = kfy;
LCC_6_to_1 = kby;
LCC_7_to_2 = kby;
LCC_8_to_3 = kby;
LCC_9_to_4 = kby;
LCC_10_to_5 = kby;

CaSSmod_1 = min(CaSS_1, 0.05); % 50 microM max on RyR rates except for k12
CaSSmod_2 = min(CaSS_2, 0.05);
CaSSmod_3 = min(CaSS_3, 0.05);
CaSSmod_4 = min(CaSS_4, 0.05);

RyR_1_to_2_CaSS_1 = k12*CaSS_1^2;       %   Set RyR channel rates (1/ms)
RyR_1_to_2_CaSS_2 = k12*CaSS_2^2;
RyR_1_to_2_CaSS_3 = k12*CaSS_3^2;
RyR_1_to_2_CaSS_4 = k12*CaSS_4^2;
RyR_2_to_1 = k21;
RyR_2_to_3_CaSS_1 = k23*CaSSmod_1^2;
RyR_2_to_3_CaSS_2 = k23*CaSSmod_2^2;
RyR_2_to_3_CaSS_3 = k23*CaSSmod_3^2;
RyR_2_to_3_CaSS_4 = k23*CaSSmod_4^2;
RyR_3_to_2_CaSS_1 = k32*k43/(k34*CaSSmod_1^2 + k43);
RyR_3_to_2_CaSS_2 = k32*k43/(k34*CaSSmod_2^2 + k43);
RyR_3_to_2_CaSS_3 = k32*k43/(k34*CaSSmod_3^2 + k43);
RyR_3_to_2_CaSS_4 = k32*k43/(k34*CaSSmod_4^2 + k43);
RyR_2_to_4_CaSS_1 = k25*CaSSmod_1^2;
RyR_2_to_4_CaSS_2 = k25*CaSSmod_2^2;
RyR_2_to_4_CaSS_3 = k25*CaSSmod_3^2;
RyR_2_to_4_CaSS_4 = k25*CaSSmod_4^2;
RyR_4_to_2_CaSS_1 = k52*k65/(k56*CaSSmod_1^2 + k65);
RyR_4_to_2_CaSS_2 = k52*k65/(k56*CaSSmod_2^2 + k65);
RyR_4_to_2_CaSS_3 = k52*k65/(k56*CaSSmod_3^2 + k65);
RyR_4_to_2_CaSS_4 = k52*k65/(k56*CaSSmod_4^2 + k65);
RyR_3_to_4_CaSS_1 = k45*k34*CaSSmod_1^2/(k34*CaSSmod_1^2 + k43);
RyR_3_to_4_CaSS_2 = k45*k34*CaSSmod_2^2/(k34*CaSSmod_2^2 + k43);
RyR_3_to_4_CaSS_3 = k45*k34*CaSSmod_3^2/(k34*CaSSmod_3^2 + k43);
RyR_3_to_4_CaSS_4 = k45*k34*CaSSmod_4^2/(k34*CaSSmod_4^2 + k43);
RyR_4_to_3_CaSS_1 = k65*k54*CaSSmod_1^2/(k56*CaSSmod_1^2 + k65);
RyR_4_to_3_CaSS_2 = k65*k54*CaSSmod_2^2/(k56*CaSSmod_2^2 + k65);
RyR_4_to_3_CaSS_3 = k65*k54*CaSSmod_3^2/(k56*CaSSmod_3^2 + k65);
RyR_4_to_3_CaSS_4 = k65*k54*CaSSmod_4^2/(k56*CaSSmod_4^2 + k65);


dy = zeros(length(y),1);    %   Evaluate state derivatives for coupled LCC-RyR model

dy( 1) = -LCC_1_to_2*y(1) +LCC_2_to_1*y(2) -LCC_1_to_4_CaSS_1*y(1) +LCC_4_to_1*y(4)  ...
         -LCC_1_to_6*y(1) +LCC_6_to_1*y(6) -RyR_1_to_2_CaSS_1*y(1) +RyR_2_to_1*y(11);
dy( 2) = +LCC_1_to_2*y(1) -LCC_2_to_1*y(2) -LCC_2_to_3*y(2) +LCC_3_to_2*y(3) ...
         -LCC_2_to_5_CaSS_1*y(2) +LCC_5_to_2*y(5) -LCC_2_to_7*y(2) +LCC_7_to_2*y(7) ...
         -RyR_1_to_2_CaSS_1*y(2) +RyR_2_to_1*y(12);
dy( 3) = +LCC_2_to_3*y(2) -LCC_3_to_2*y(3) -LCC_3_to_8*y(3) +LCC_8_to_3*y(8) ...
         -RyR_1_to_2_CaSS_2*y(3) +RyR_2_to_1*y(13);
dy( 4) = +LCC_1_to_4_CaSS_1*y(1) -LCC_4_to_1*y(4) -LCC_4_to_5*y(4) +LCC_5_to_4*y(5) ...
         -LCC_4_to_9*y(4) +LCC_9_to_4*y(9) -RyR_1_to_2_CaSS_1*y(4) +RyR_2_to_1*y(14);
dy( 5) = +LCC_2_to_5_CaSS_1*y(2) -LCC_5_to_2*y(5) +LCC_4_to_5*y(4) -LCC_5_to_4*y(5) ...
         -LCC_5_to_10*y(5) +LCC_10_to_5*y(10) -RyR_1_to_2_CaSS_1*y(5) +RyR_2_to_1*y(15);
dy( 6) = +LCC_1_to_6*y(1) -LCC_6_to_1*y(6) -LCC_6_to_7*y(6) +LCC_7_to_6*y(7) ...
         -LCC_6_to_9_CaSS_1*y(6) +LCC_9_to_6*y(9) -RyR_1_to_2_CaSS_1*y(6) +RyR_2_to_1*y(16);
dy( 7) = +LCC_2_to_7*y(2) -LCC_7_to_2*y(7) +LCC_6_to_7*y(6) -LCC_7_to_6*y(7) ...
         -LCC_7_to_8*y(7) +LCC_8_to_7*y(8) -LCC_7_to_10_CaSS_1*y(7) +LCC_10_to_7*y(10) ...
         -RyR_1_to_2_CaSS_1*y(7) +RyR_2_to_1*y(17);
dy( 8) = +LCC_3_to_8*y(3) -LCC_8_to_3*y(8) +LCC_7_to_8*y(7) -LCC_8_to_7*y(8) ...
         -RyR_1_to_2_CaSS_1*y(8) +RyR_2_to_1*y(18);
dy( 9) = +LCC_4_to_9*y(4) -LCC_9_to_4*y(9) +LCC_6_to_9_CaSS_1*y(6) -LCC_9_to_6*y(9) ...
         -LCC_9_to_10*y(9) +LCC_10_to_9*y(10) -RyR_1_to_2_CaSS_1*y(9) +RyR_2_to_1*y(19);
dy(10) = +LCC_5_to_10*y(5) -LCC_10_to_5*y(10) +LCC_7_to_10_CaSS_1*y(7) -LCC_10_to_7*y(10) ...
         +LCC_9_to_10*y(9) -LCC_10_to_9*y(10) -RyR_1_to_2_CaSS_1*y(10) +RyR_2_to_1*y(20);
dy(11) = +RyR_1_to_2_CaSS_1*y(1) -RyR_2_to_1*y(11) -LCC_1_to_2*y(11) +LCC_2_to_1*y(12) ...
         -LCC_1_to_4_CaSS_1*y(11) +LCC_4_to_1*y(14) -LCC_1_to_6*y(11) +LCC_6_to_1*y(16) ...
         -RyR_2_to_3_CaSS_1*y(11) +RyR_3_to_2_CaSS_3*y(21) -RyR_2_to_4_CaSS_1*y(11) +RyR_4_to_2_CaSS_1*y(31);
dy(12) = +RyR_1_to_2_CaSS_1*y(2) -RyR_2_to_1*y(12) +LCC_1_to_2*y(11) -LCC_2_to_1*y(12) ...
         -LCC_2_to_3*y(12) +LCC_3_to_2*y(13) -LCC_2_to_5_CaSS_1*y(12) +LCC_5_to_2*y(15) ...
         -LCC_2_to_7*y(12) +LCC_7_to_2*y(17) -RyR_2_to_3_CaSS_1*y(12) +RyR_3_to_2_CaSS_3*y(22) ...
         -RyR_2_to_4_CaSS_1*y(12) +RyR_4_to_2_CaSS_1*y(32);
dy(13) = +RyR_1_to_2_CaSS_2*y(3) -RyR_2_to_1*y(13) +LCC_2_to_3*y(12) -LCC_3_to_2*y(13) ...
         -LCC_3_to_8*y(13) +LCC_8_to_3*y(18) -RyR_2_to_3_CaSS_2*y(13) +RyR_3_to_2_CaSS_4*y(23) ...
         -RyR_2_to_4_CaSS_2*y(13) +RyR_4_to_2_CaSS_2*y(33);
dy(14) = +RyR_1_to_2_CaSS_1*y(4) -RyR_2_to_1*y(14) +LCC_1_to_4_CaSS_1*y(11) -LCC_4_to_1*y(14) ...
         -LCC_4_to_5*y(14) +LCC_5_to_4*y(15) -LCC_4_to_9*y(14) +LCC_9_to_4*y(19) ...
         -RyR_2_to_3_CaSS_1*y(14) +RyR_3_to_2_CaSS_3*y(24) -RyR_2_to_4_CaSS_1*y(14) +RyR_4_to_2_CaSS_1*y(34);
dy(15) = +RyR_1_to_2_CaSS_1*y(5) -RyR_2_to_1*y(15) +LCC_2_to_5_CaSS_1*y(12) -LCC_5_to_2*y(15) ...
         +LCC_4_to_5*y(14) -LCC_5_to_4*y(15) -LCC_5_to_10*y(15) +LCC_10_to_5*y(20) ...
         -RyR_2_to_3_CaSS_1*y(15) +RyR_3_to_2_CaSS_3*y(25) -RyR_2_to_4_CaSS_1*y(15) +RyR_4_to_2_CaSS_1*y(35);
dy(16) = +RyR_1_to_2_CaSS_1*y(6) -RyR_2_to_1*y(16) +LCC_1_to_6*y(11) -LCC_6_to_1*y(16) ...
         -LCC_6_to_7*y(16) +LCC_7_to_6*y(17) -LCC_6_to_9_CaSS_1*y(16) +LCC_9_to_6*y(19) ...
         -RyR_2_to_3_CaSS_1*y(16) +RyR_3_to_2_CaSS_3*y(26) -RyR_2_to_4_CaSS_1*y(16) +RyR_4_to_2_CaSS_1*y(36);
dy(17) = +RyR_1_to_2_CaSS_1*y(7) -RyR_2_to_1*y(17) +LCC_2_to_7*y(12) -LCC_7_to_2*y(17) ...
         +LCC_6_to_7*y(16) -LCC_7_to_6*y(17) -LCC_7_to_8*y(17) +LCC_8_to_7*y(18) ...
         -LCC_7_to_10_CaSS_1*y(17) +LCC_10_to_7*y(20) -RyR_2_to_3_CaSS_1*y(17) +RyR_3_to_2_CaSS_3*y(27) ...
         -RyR_2_to_4_CaSS_1*y(17) +RyR_4_to_2_CaSS_1*y(37);
dy(18) = +RyR_1_to_2_CaSS_1*y(8) -RyR_2_to_1*y(18) +LCC_3_to_8*y(13) -LCC_8_to_3*y(18) ...
         +LCC_7_to_8*y(17) -LCC_8_to_7*y(18) -RyR_2_to_3_CaSS_1*y(18) +RyR_3_to_2_CaSS_3*y(28) ...
         -RyR_2_to_4_CaSS_1*y(18) +RyR_4_to_2_CaSS_1*y(38);
dy(19) = +RyR_1_to_2_CaSS_1*y(9) -RyR_2_to_1*y(19) +LCC_4_to_9*y(14) -LCC_9_to_4*y(19) ...
         +LCC_6_to_9_CaSS_1*y(16) -LCC_9_to_6*y(19) -LCC_9_to_10*y(19) +LCC_10_to_9*y(20) ...
         -RyR_2_to_3_CaSS_1*y(19) +RyR_3_to_2_CaSS_3*y(29) -RyR_2_to_4_CaSS_1*y(19) +RyR_4_to_2_CaSS_1*y(39);
dy(20) = +RyR_1_to_2_CaSS_1*y(10) -RyR_2_to_1*y(20) +LCC_5_to_10*y(15) -LCC_10_to_5*y(20) ...
         +LCC_7_to_10_CaSS_1*y(17) -LCC_10_to_7*y(20) +LCC_9_to_10*y(19) -LCC_10_to_9*y(20) ...
         -RyR_2_to_3_CaSS_1*y(20) +RyR_3_to_2_CaSS_3*y(30) -RyR_2_to_4_CaSS_1*y(20) +RyR_4_to_2_CaSS_1*y(40);
dy(21) = +RyR_2_to_3_CaSS_1*y(11) -RyR_3_to_2_CaSS_3*y(21) -LCC_1_to_2*y(21) +LCC_2_to_1*y(22) ...
         -LCC_1_to_4_CaSS_3*y(21) +LCC_4_to_1*y(24) -LCC_1_to_6*y(21) +LCC_6_to_1*y(26) ...
         -RyR_3_to_4_CaSS_3*y(21) +RyR_4_to_3_CaSS_1*y(31);
dy(22) = +RyR_2_to_3_CaSS_1*y(12) -RyR_3_to_2_CaSS_3*y(22) +LCC_1_to_2*y(21) -LCC_2_to_1*y(22) ...
         -LCC_2_to_3*y(22) +LCC_3_to_2*y(23) -LCC_2_to_5_CaSS_3*y(22) +LCC_5_to_2*y(25) ...
         -LCC_2_to_7*y(22) +LCC_7_to_2*y(27) -RyR_3_to_4_CaSS_3*y(22) +RyR_4_to_3_CaSS_1*y(32);
dy(23) = +RyR_2_to_3_CaSS_2*y(13) -RyR_3_to_2_CaSS_4*y(23) +LCC_2_to_3*y(22) -LCC_3_to_2*y(23) ...
         -LCC_3_to_8*y(23) +LCC_8_to_3*y(28) -RyR_3_to_4_CaSS_4*y(23) +RyR_4_to_3_CaSS_2*y(33);
dy(24) = +RyR_2_to_3_CaSS_1*y(14) -RyR_3_to_2_CaSS_3*y(24) +LCC_1_to_4_CaSS_3*y(21) -LCC_4_to_1*y(24) ...
         -LCC_4_to_5*y(24) +LCC_5_to_4*y(25) -LCC_4_to_9*y(24) +LCC_9_to_4*y(29) ...
         -RyR_3_to_4_CaSS_3*y(24) +RyR_4_to_3_CaSS_1*y(34);
dy(25) = +RyR_2_to_3_CaSS_1*y(15) -RyR_3_to_2_CaSS_3*y(25) +LCC_2_to_5_CaSS_3*y(22) -LCC_5_to_2*y(25) ...
         +LCC_4_to_5*y(24) -LCC_5_to_4*y(25) -LCC_5_to_10*y(25) +LCC_10_to_5*y(30) ...
         -RyR_3_to_4_CaSS_3*y(25) +RyR_4_to_3_CaSS_1*y(35);
dy(26) = +RyR_2_to_3_CaSS_1*y(16) -RyR_3_to_2_CaSS_3*y(26) +LCC_1_to_6*y(21) -LCC_6_to_1*y(26) ...
         -LCC_6_to_7*y(26) +LCC_7_to_6*y(27) -LCC_6_to_9_CaSS_3*y(26) +LCC_9_to_6*y(29) ...
         -RyR_3_to_4_CaSS_3*y(26) +RyR_4_to_3_CaSS_1*y(36);
dy(27) = +RyR_2_to_3_CaSS_1*y(17) -RyR_3_to_2_CaSS_3*y(27) +LCC_2_to_7*y(22) -LCC_7_to_2*y(27) ...
         +LCC_6_to_7*y(26) -LCC_7_to_6*y(27) -LCC_7_to_8*y(27) +LCC_8_to_7*y(28) ...
         -LCC_7_to_10_CaSS_3*y(27) +LCC_10_to_7*y(30) -RyR_3_to_4_CaSS_3*y(27) +RyR_4_to_3_CaSS_1*y(37);
dy(28) = +RyR_2_to_3_CaSS_1*y(18) -RyR_3_to_2_CaSS_3*y(28) +LCC_3_to_8*y(23) -LCC_8_to_3*y(28) ...
         +LCC_7_to_8*y(27) -LCC_8_to_7*y(28) -RyR_3_to_4_CaSS_3*y(28) +RyR_4_to_3_CaSS_1*y(38);
dy(29) = +RyR_2_to_3_CaSS_1*y(19) -RyR_3_to_2_CaSS_3*y(29) +LCC_4_to_9*y(24) -LCC_9_to_4*y(29) ...
         +LCC_6_to_9_CaSS_3*y(26) -LCC_9_to_6*y(29) -LCC_9_to_10*y(29) +LCC_10_to_9*y(30) ...
         -RyR_3_to_4_CaSS_3*y(29) +RyR_4_to_3_CaSS_1*y(39);
dy(30) = +RyR_2_to_3_CaSS_1*y(20) -RyR_3_to_2_CaSS_3*y(30) +LCC_5_to_10*y(25) -LCC_10_to_5*y(30) ...
         +LCC_7_to_10_CaSS_3*y(27) -LCC_10_to_7*y(30) +LCC_9_to_10*y(29) -LCC_10_to_9*y(30) ...
         -RyR_3_to_4_CaSS_3*y(30) +RyR_4_to_3_CaSS_1*y(40);
dy(31) = +RyR_2_to_4_CaSS_1*y(11) -RyR_4_to_2_CaSS_1*y(31) +RyR_3_to_4_CaSS_3*y(21) -RyR_4_to_3_CaSS_1*y(31) ...
         -LCC_1_to_2*y(31) +LCC_2_to_1*y(32) -LCC_1_to_4_CaSS_1*y(31) +LCC_4_to_1*y(34) ...
         -LCC_1_to_6*y(31) +LCC_6_to_1*y(36);
dy(32) = +RyR_2_to_4_CaSS_1*y(12) -RyR_4_to_2_CaSS_1*y(32) +RyR_3_to_4_CaSS_3*y(22) -RyR_4_to_3_CaSS_1*y(32) ...
         +LCC_1_to_2*y(31) -LCC_2_to_1*y(32) -LCC_2_to_3*y(32) +LCC_3_to_2*y(33) ...
         -LCC_2_to_5_CaSS_1*y(32) +LCC_5_to_2*y(35) -LCC_2_to_7*y(32) +LCC_7_to_2*y(37);
dy(33) = +RyR_2_to_4_CaSS_2*y(13) -RyR_4_to_2_CaSS_2*y(33) +RyR_3_to_4_CaSS_4*y(23) -RyR_4_to_3_CaSS_2*y(33) ...
         +LCC_2_to_3*y(32) -LCC_3_to_2*y(33) -LCC_3_to_8*y(33) +LCC_8_to_3*y(38);
dy(34) = +RyR_2_to_4_CaSS_1*y(14) -RyR_4_to_2_CaSS_1*y(34) +RyR_3_to_4_CaSS_3*y(24) -RyR_4_to_3_CaSS_1*y(34) ...
         +LCC_1_to_4_CaSS_1*y(31) -LCC_4_to_1*y(34) -LCC_4_to_5*y(34) +LCC_5_to_4*y(35) ...
         -LCC_4_to_9*y(34) +LCC_9_to_4*y(39);
dy(35) = +RyR_2_to_4_CaSS_1*y(15) -RyR_4_to_2_CaSS_1*y(35) +RyR_3_to_4_CaSS_3*y(25) -RyR_4_to_3_CaSS_1*y(35) ...
         +LCC_2_to_5_CaSS_1*y(32) -LCC_5_to_2*y(35) +LCC_4_to_5*y(34) -LCC_5_to_4*y(35) ...
         -LCC_5_to_10*y(35) +LCC_10_to_5*y(40);
dy(36) = +RyR_2_to_4_CaSS_1*y(16) -RyR_4_to_2_CaSS_1*y(36) +RyR_3_to_4_CaSS_3*y(26) -RyR_4_to_3_CaSS_1*y(36) ...
         +LCC_1_to_6*y(31) -LCC_6_to_1*y(36) -LCC_6_to_7*y(36) +LCC_7_to_6*y(37) ...
         -LCC_6_to_9_CaSS_1*y(36) +LCC_9_to_6*y(39);
dy(37) = +RyR_2_to_4_CaSS_1*y(17) -RyR_4_to_2_CaSS_1*y(37) +RyR_3_to_4_CaSS_3*y(27) -RyR_4_to_3_CaSS_1*y(37) ...
         +LCC_2_to_7*y(32) -LCC_7_to_2*y(37) +LCC_6_to_7*y(36) -LCC_7_to_6*y(37) ...
         -LCC_7_to_8*y(37) +LCC_8_to_7*y(38) -LCC_7_to_10_CaSS_1*y(37) +LCC_10_to_7*y(40);
dy(38) = +RyR_2_to_4_CaSS_1*y(18) -RyR_4_to_2_CaSS_1*y(38) +RyR_3_to_4_CaSS_3*y(28) -RyR_4_to_3_CaSS_1*y(38) ...
         +LCC_3_to_8*y(33) -LCC_8_to_3*y(38) +LCC_7_to_8*y(37) -LCC_8_to_7*y(38);
dy(39) = +RyR_2_to_4_CaSS_1*y(19) -RyR_4_to_2_CaSS_1*y(39) +RyR_3_to_4_CaSS_3*y(29) -RyR_4_to_3_CaSS_1*y(39) ...
         +LCC_4_to_9*y(34) -LCC_9_to_4*y(39) +LCC_6_to_9_CaSS_1*y(36) -LCC_9_to_6*y(39) ...
         -LCC_9_to_10*y(39) +LCC_10_to_9*y(40);
dy(40) = +RyR_2_to_4_CaSS_1*y(20) -RyR_4_to_2_CaSS_1*y(40) +RyR_3_to_4_CaSS_3*y(30) -RyR_4_to_3_CaSS_1*y(40) ...
         +LCC_5_to_10*y(35) -LCC_10_to_5*y(40) +LCC_7_to_10_CaSS_1*y(37) -LCC_10_to_7*y(40) ...
         +LCC_9_to_10*y(39) -LCC_10_to_9*y(40);


PCaSS_2 = y(3)+y(13)+y(33);         %   Evaluate occupancy probability for each open-closed configuartion
PCaSS_3 = y(21)+y(22)+y(24)+y(25)+y(26)+y(27)+y(28)+y(29)+y(30);
PCaSS_4 = y(23);
PCaSS_1 = 1 - PCaSS_2 - PCaSS_3 - PCaSS_4;

                                    %   Evaluate fluxes/currents across LCC-RyR boundary
ICaL = NCaRU/Csa*PCa*4*VFsq_over_RT/(exp(2*VF_over_RT)-1) ...
    * ( PCaSS_2*(CaSS_2*exp(2*VF_over_RT)-0.341*Cao) + PCaSS_4*(CaSS_4*exp(2*VF_over_RT)-0.341*Cao) );

JCaL = -ICaL*Csa/(2*Faraday*1000*VSS);
Ito2 = Ito2scale*NCaRU/Csa*PCl*VFsq_over_RT*(Cli*exp(-VF_over_RT)-Clo)/(exp(-VF_over_RT) - 1) ...
        * ( PCaSS_4/(1 + KdIto2/CaSS_4) + PCaSS_3/(1 + KdIto2/CaSS_3) ...
            + PCaSS_2/(1 + KdIto2/CaSS_2) + PCaSS_1/(1 + KdIto2/CaSS_1) );

JRyR = NCaRU*JRyRmax * ( PCaSS_4*(CaSR - CaSS_4) + PCaSS_3*(CaSR - CaSS_3) );
Jss2b = NCaRU*r_xfer * ( PCaSS_4*(CaSS_4 - Cai) + PCaSS_3*(CaSS_3 - Cai) ...
                            + PCaSS_2*(CaSS_2 - Cai) + PCaSS_1*(CaSS_1 - Cai) );
CaSSavg = PCaSS_1*CaSS_1 + PCaSS_2*CaSS_2 + PCaSS_3*CaSS_3 + PCaSS_4*CaSS_4;        % Average dyadic Ca concentration


                                    %   Evaluate currents and state derivatives of global model 
ENa = RT_over_F*log(Nao/Nai);
EK =  RT_over_F*log(Ko/Ki);
EKs = RT_over_F*log((Ko+0.01833*Nao)/(Ki+0.01833*Nai));
ECa = 0.5*RT_over_F*log(Cao/max(1e-10,Cai));

%%%% Markov model for WT INa including INaL - Clancy 2003%%%%%%
a11= 3.802/(0.1027*exp(-V/17.0)+0.20*exp(-V/150));
a12=  ( 3.802/(0.1027*exp(-V/15.0)+0.23*exp(-V/150)));
a13=3.802/(0.1027*exp(-V/12.0)+0.250*exp(-V/150));
b11= 0.4*exp(-V/20.3); %Clancy et al, 2003 (WT/I1768V)
b12= 0.4*exp(-(V-5)/20.3); %Clancy et al, 2003 (WT/I1768V)
b13= 0.4*exp(-(V-10)/20.3)/4.5; %Clancy et al, 2003 (WT/I1768V)
a3= (3.7933e-7*exp(-V/7.7))*5*i768v;%Clancy et al, 2003 (WT/I1768V)
b3= (0.0084+.00002*V);
a2= ((9.178*exp(V/29.68))/4.5);%Clancy et al, 2003 (WT/I1768V)
b2= ((a13*a2*a3)/(b13*b3));	
a4= (a2/100)*1.5; %Clancy et al, 2003 (WT/I1768V)
b4 = a3*5*i768v;%Clancy et al, 2003 (WT/I1768V)
a5= (a2/95000)*1.5; %Clancy et al, 2003 (WT/I1768V)
b5 = (a3/50)*5*i768v;%Clancy et al, 2003 (WT/I1768V)

mu1= 1.0e-7;
mu2= 3.0e-4;

%Upper states
dP_UO = (P_UC1*a13+P_UIF*b2+P_LO*mu2)- (P_UO*(b13+a2+mu1));
dP_UIF = (P_IC2*a12+P_UO*a2+P_UC1*b3+P_UIM1*b4)- (P_UIF*(b2+a3+a4+b12));
dP_UC2 = (P_IC2*a3+P_UC1*b12+P_UC3*a11+P_LC2*mu2)-(P_UC2*(a12+b11+mu1+b3));
dP_UC3 = P_IC3*a3+P_UC2*b11+P_LC3*mu2-(P_UC3*(a11+mu1+b3));
dP_UIM1 = (P_UIF*a4+P_UIM2*b5)-P_UIM1*(b4+a5);
dP_UC1 = P_UC2*a12+P_UO*b13+P_LC1*mu2+P_UIF*a3-(P_UC1*(a13+b12+mu1+b3));
dP_IC3 = (P_UC3*b3+P_IC2*b11)-P_IC3*(a11+a3);
dP_IC2 = (P_UC2*b3+P_IC3*a11+P_UIF*b12)-P_IC2*(a12+b11+a3);
dP_UIM2 = P_UIM1*a5-P_UIM2*b5;

%Lower states
dP_LC3 = P_UC3*mu1+P_LC2*b11-(P_LC3*(a11+mu2));
dP_LC2 = P_UC2*mu1+P_LC3*a11+P_LC1*b12-(P_LC2*(a12+mu2+b11));
dP_LC1 = P_UC1*mu1+P_LC2*a12+P_LO*b13-(P_LC1*(a13+mu2+b12));
dP_LO = (P_LC1*a13+P_UO*mu1)-(P_LO*(b13+mu2));


GNa = 4.6; % mS/uF
INa = GNa*(P_UO)*(V-ENa); % INa 

INal = GNal*(P_LO)*(V-ENa); % INal

fKo = (Ko/4)^0.5;
IKr = GKr*fKo*OHerg*(V-EK);

IKs = GKs*(xKs^2)*(V-EKs);

IKv43 = GKv43*(OKv43)*(V-EK);

a1 =  Ki*exp(VF_over_RT)-Ko;
a2 =  exp(VF_over_RT)-1;
IKv14_K = PKv14*OKv14*VFsq_over_RT*(a1/a2);
a1 =  Nai*exp(VF_over_RT)-Nao;
IKv14_Na = 0.02*PKv14*OKv14*VFsq_over_RT*(a1/a2);
IKv14 = IKv14_K + IKv14_Na;
Ito1 = IKv43 + IKv14;

K1_inf = 1/(0.94+exp(1.76/RT_over_F*(V-EK)));
IK1 = GK1*Ko/(Ko+KmK1)*K1_inf*(V-EK);

INab = GNab*(V-ENa);

KpV = 1/(1+exp((7.488-V)/5.98));
IKp = GKp*KpV*(V-EK);				
      
sigma = (exp(Nao/67.3)-1)/7;
a1 = 1+0.1245*exp(-0.1*VF_over_RT);
a2 = 0.0365*sigma*exp(-VF_over_RT);
fNaK = 1/(a1+a2);
a1 = Ko/(Ko+KmKo);
a2 = 1+(KmNai/Nai)^2;
INaK = INaKmax*fNaK*(a1/a2);


a1 = exp(eta*VF_over_RT)*Nai^3*Cao;
a2 = exp((eta-1)*VF_over_RT)*Nao^3*Cai;
a3 = 1+ksat*exp((eta-1)*VF_over_RT);
a4 = KmCa+Cao;
a5 = 5000/(KmNa^3+Nao^3);
INaCa = kNaCa*a5*(a1-a2)/(a4*a3);

ICab = GCab*(V-ECa);

IpCa = IpCamax*Cai/(KmpCa+Cai);


C1H_to_C2H = T_Const_HERG*A0_HERG*exp(B0_HERG*V);
C2H_to_C1H = T_Const_HERG*A1_HERG*exp(B1_HERG*V);
C3H_to_OH =  T_Const_HERG*A2_HERG*exp(B2_HERG*V);
OH_to_C3H =  T_Const_HERG*A3_HERG*exp(B3_HERG*V);
OH_to_IH =   T_Const_HERG*A4_HERG*exp(B4_HERG*V);
IH_to_OH =   T_Const_HERG*A5_HERG*exp(B5_HERG*V);
C3H_to_IH =  T_Const_HERG*A6_HERG*exp(B6_HERG*V);
IH_to_C3H =  (OH_to_C3H*IH_to_OH*C3H_to_IH)/(C3H_to_OH*OH_to_IH);
dC1Herg = C2H_to_C1H * C2Herg - C1H_to_C2H * C1Herg;
a1 = C1H_to_C2H * C1Herg + C3H_to_C2H * C3Herg;
a2 = (C2H_to_C1H + C2H_to_C3H) * C2Herg;
dC2Herg = a1-a2;
a1 = C2H_to_C3H*C2Herg + OH_to_C3H*OHerg + IH_to_C3H*IHerg;
a2 = (C3H_to_IH + C3H_to_OH + C3H_to_C2H) * C3Herg;
dC3Herg = a1-a2;
a1 = C3H_to_OH * C3Herg + IH_to_OH * IHerg;
a2 = (OH_to_C3H + OH_to_IH) * OHerg;
dOHerg = a1-a2;
a1 = C3H_to_IH * C3Herg + OH_to_IH * OHerg;
a2 = (IH_to_C3H + IH_to_OH) * IHerg;
dIHerg = a1-a2;

xKs_inf = 1/(1+exp(-(V-24.7)/13.6));
a1 = 7.19e-5*(V-10)/(1-exp(-0.148*(V-10)));
a2 = 1.31e-4*(V-10)/(exp(0.0687*(V-10))-1);
tau_xKs = 1/(a1+a2);
dxKs = (xKs_inf-xKs)/tau_xKs;
	
fb = (Cai/Kfb)^Nfb;
rb = (CaSR/Krb)^Nrb;
Jup = (vmaxf*fb - vmaxr*rb)/(1 + fb + rb);
 
dLTRPNCa = kltrpn_plus*Cai*(1 - LTRPNCa) - kltrpn_minus * LTRPNCa;
dHTRPNCa = khtrpn_plus*Cai*(1 - HTRPNCa) - khtrpn_minus * HTRPNCa;
Jtrpn = LTRPNtot*dLTRPNCa+HTRPNtot*dHTRPNCa;

beta_i = 1/(1+CMDNtot*KmCMDN/(Cai+KmCMDN)^2+EGTAtot*KmEGTA/(Cai+KmEGTA)^2);
beta_SR = 1/(1+CSQNtot*KmCSQN/(CaSR+KmCSQN)^2);

a1 = Csa/(Vmyo*Faraday*1000);
dNai = -( INal + INa+INab+3*(INaCa+INaK)+ IKv14_Na)*a1;

a3 = IKr+IKs+IK1+IKp;
dKi = -( a3-2*INaK+Ito1 )*a1;

a3 = ICab-2*INaCa+IpCa;
dCai = beta_i*(Jss2b*VSS/Vmyo-Jup-Jtrpn - a3*0.5*a1);

dCaSR = beta_SR*(Jup*Vmyo/VSR - JRyR*VSS/VSR);


a1 = INal+ INa+ICaL+IKr+IKs;
a2 = IK1+IKp+INaCa+INaK+Ito1+Ito2;
a3 = IpCa+ICab+INab;
Itot = a1+a2+a3+Istim;
dV = -Itot;

alpha_act43 = alphaa0Kv43*exp(aaKv43*V);
beta_act43  = betaa0Kv43*exp(-baKv43*V);

alpha_inact43 = alphai0Kv43*exp(-aiKv43*V);
beta_inact43  = betai0Kv43*exp(biKv43*V);
alpha_act43_Ca = k_alphaoi*Cai;
beta_act43_Ca = k_betaoi*exp(-betaoi*V);


C0Kv43_to_C1Kv43 = 4*alpha_act43;
C1Kv43_to_C2Kv43 = 3*alpha_act43;
C2Kv43_to_C3Kv43 = 2*alpha_act43;
C3Kv43_to_OKv43  =   alpha_act43;

CI0Kv43_to_CI1Kv43 = 4*b1Kv43*alpha_act43;
CI1Kv43_to_CI2Kv43 = 3*b2Kv43*alpha_act43/b1Kv43;
CI2Kv43_to_CI3Kv43 = 2*b3Kv43*alpha_act43/b2Kv43;
CI3Kv43_to_OIKv43  =   b4Kv43*alpha_act43/b3Kv43;

C1Kv43_to_C0Kv43 =   beta_act43;
C2Kv43_to_C1Kv43 = 2*beta_act43;
C3Kv43_to_C2Kv43 = 3*beta_act43;
OKv43_to_C3Kv43  = 4*beta_act43;

CI1Kv43_to_CI0Kv43 =          beta_act43/f1Kv43;
CI2Kv43_to_CI1Kv43 = 2*f1Kv43*beta_act43/f2Kv43;
CI3Kv43_to_CI2Kv43 = 3*f2Kv43*beta_act43/f3Kv43;
OIKv43_to_CI3Kv43  = 4*f3Kv43*beta_act43/f4Kv43;

C0Kv43_to_CI0Kv43 = beta_inact43;
C1Kv43_to_CI1Kv43 = f1Kv43*beta_inact43;
C2Kv43_to_CI2Kv43 = f2Kv43*beta_inact43;
C3Kv43_to_CI3Kv43 = f3Kv43*beta_inact43;
OKv43_to_OIKv43   = f4Kv43*beta_inact43;

OKv43_to_OIKv43Ca   = alpha_act43_Ca;
OIKv43Ca_to_OKv43   = beta_act43_Ca;

CI0Kv43_to_C0Kv43 = alpha_inact43;
CI1Kv43_to_C1Kv43 = alpha_inact43/b1Kv43;
CI2Kv43_to_C2Kv43 = alpha_inact43/b2Kv43;
CI3Kv43_to_C3Kv43 = alpha_inact43/b3Kv43;
OIKv43_to_OKv43   = alpha_inact43/b4Kv43;

a1 = (C0Kv43_to_C1Kv43+C0Kv43_to_CI0Kv43)*C0Kv43;
a2 = C1Kv43_to_C0Kv43*C1Kv43 + CI0Kv43_to_C0Kv43*CI0Kv43;
dC0Kv43 = a2 - a1;

a1 = (C1Kv43_to_C2Kv43+C1Kv43_to_C0Kv43+C1Kv43_to_CI1Kv43)*C1Kv43;
a2 = C2Kv43_to_C1Kv43*C2Kv43 + CI1Kv43_to_C1Kv43*CI1Kv43 + C0Kv43_to_C1Kv43*C0Kv43;
dC1Kv43 = a2 - a1;

a1 = (C2Kv43_to_C3Kv43+C2Kv43_to_C1Kv43+C2Kv43_to_CI2Kv43)*C2Kv43;
a2 = C3Kv43_to_C2Kv43*C3Kv43 + CI2Kv43_to_C2Kv43*CI2Kv43 + C1Kv43_to_C2Kv43*C1Kv43;
dC2Kv43 = a2 - a1;

a1 = (C3Kv43_to_OKv43+C3Kv43_to_C2Kv43+C3Kv43_to_CI3Kv43)*C3Kv43;
a2 = OKv43_to_C3Kv43*OKv43 + CI3Kv43_to_C3Kv43*CI3Kv43 + C2Kv43_to_C3Kv43*C2Kv43;
dC3Kv43 = a2 - a1;

a1 = (OKv43_to_C3Kv43+OKv43_to_OIKv43+ OKv43_to_OIKv43Ca )*OKv43;
a2 = C3Kv43_to_OKv43*C3Kv43 + OIKv43_to_OKv43*OIKv43 + OIKv43Ca_to_OKv43*OIKv43Ca;
dOKv43 = a2 - a1;

a1 = (CI0Kv43_to_C0Kv43+CI0Kv43_to_CI1Kv43)*CI0Kv43;
a2 = C0Kv43_to_CI0Kv43*C0Kv43 + CI1Kv43_to_CI0Kv43*CI1Kv43;
dCI0Kv43 = a2 - a1;

a1 = (CI1Kv43_to_CI2Kv43+CI1Kv43_to_C1Kv43+CI1Kv43_to_CI0Kv43)*CI1Kv43;
a2 = CI2Kv43_to_CI1Kv43*CI2Kv43 + C1Kv43_to_CI1Kv43*C1Kv43 + CI0Kv43_to_CI1Kv43*CI0Kv43;
dCI1Kv43 = a2 - a1;

a1 = (CI2Kv43_to_CI3Kv43+CI2Kv43_to_C2Kv43+CI2Kv43_to_CI1Kv43)*CI2Kv43;
a2 = CI3Kv43_to_CI2Kv43*CI3Kv43 + C2Kv43_to_CI2Kv43*C2Kv43 + CI1Kv43_to_CI2Kv43*CI1Kv43;
dCI2Kv43 = a2 - a1;

a1 = (CI3Kv43_to_OIKv43+CI3Kv43_to_C3Kv43+CI3Kv43_to_CI2Kv43)*CI3Kv43;
a2 = OIKv43_to_CI3Kv43*OIKv43 + C3Kv43_to_CI3Kv43*C3Kv43 + CI2Kv43_to_CI3Kv43*CI2Kv43;
dCI3Kv43 = a2 - a1;

a1 = (OIKv43_to_OKv43+OIKv43_to_CI3Kv43)*OIKv43;
a2 = OKv43_to_OIKv43*OKv43 + CI3Kv43_to_OIKv43*CI3Kv43;
dOIKv43 = a2 - a1;



a1 = (OIKv43Ca_to_OKv43)*OIKv43Ca;
a2 = OKv43_to_OIKv43Ca*OKv43;
dOIKv43Ca = a2-a1 ;

alpha_act14 = alphaa0Kv14*exp(aaKv14*V);
beta_act14  = betaa0Kv14*exp(-baKv14*V);
alpha_inact14 = alphai0Kv14;
beta_inact14  = betai0Kv14;

C0Kv14_to_C1Kv14 = 4*alpha_act14;
C1Kv14_to_C2Kv14 = 3*alpha_act14;
C2Kv14_to_C3Kv14 = 2*alpha_act14;
C3Kv14_to_OKv14  =   alpha_act14;

CI0Kv14_to_CI1Kv14 = 4*b1Kv14*alpha_act14;
CI1Kv14_to_CI2Kv14 = 3*b2Kv14*alpha_act14/b1Kv14;
CI2Kv14_to_CI3Kv14 = 2*b3Kv14*alpha_act14/b2Kv14;
CI3Kv14_to_OIKv14  =   b4Kv14*alpha_act14/b3Kv14;

C1Kv14_to_C0Kv14 =   beta_act14;
C2Kv14_to_C1Kv14 = 2*beta_act14;
C3Kv14_to_C2Kv14 = 3*beta_act14;
OKv14_to_C3Kv14  = 4*beta_act14;

CI1Kv14_to_CI0Kv14 =          beta_act14/f1Kv14;
CI2Kv14_to_CI1Kv14 = 2*f1Kv14*beta_act14/f2Kv14;
CI3Kv14_to_CI2Kv14 = 3*f2Kv14*beta_act14/f3Kv14;
OIKv14_to_CI3Kv14  = 4*f3Kv14*beta_act14/f4Kv14;

C0Kv14_to_CI0Kv14 = beta_inact14;
C1Kv14_to_CI1Kv14 = f1Kv14*beta_inact14;
C2Kv14_to_CI2Kv14 = f2Kv14*beta_inact14;
C3Kv14_to_CI3Kv14 = f3Kv14*beta_inact14;
OKv14_to_OIKv14   = f4Kv14*beta_inact14;

CI0Kv14_to_C0Kv14 = alpha_inact14;
CI1Kv14_to_C1Kv14 = alpha_inact14/b1Kv14;
CI2Kv14_to_C2Kv14 = alpha_inact14/b2Kv14;
CI3Kv14_to_C3Kv14 = alpha_inact14/b3Kv14;
OIKv14_to_OKv14   = alpha_inact14/b4Kv14;

a1 = (C0Kv14_to_C1Kv14+C0Kv14_to_CI0Kv14)*C0Kv14;
a2 = C1Kv14_to_C0Kv14*C1Kv14 + CI0Kv14_to_C0Kv14*CI0Kv14;
dC0Kv14 = a2 - a1;

a1 = (C1Kv14_to_C2Kv14+C1Kv14_to_C0Kv14+C1Kv14_to_CI1Kv14)*C1Kv14;
a2 = C2Kv14_to_C1Kv14*C2Kv14 + CI1Kv14_to_C1Kv14*CI1Kv14 + C0Kv14_to_C1Kv14*C0Kv14;
dC1Kv14 = a2 - a1;

a1 = (C2Kv14_to_C3Kv14+C2Kv14_to_C1Kv14+C2Kv14_to_CI2Kv14)*C2Kv14;
a2 = C3Kv14_to_C2Kv14*C3Kv14 + CI2Kv14_to_C2Kv14*CI2Kv14 + C1Kv14_to_C2Kv14*C1Kv14;
dC2Kv14 = a2 - a1;

a1 = (C3Kv14_to_OKv14+C3Kv14_to_C2Kv14+C3Kv14_to_CI3Kv14)*C3Kv14;
a2 = OKv14_to_C3Kv14*OKv14 + CI3Kv14_to_C3Kv14*CI3Kv14 + C2Kv14_to_C3Kv14*C2Kv14;
dC3Kv14 = a2 - a1;

a1 = (OKv14_to_C3Kv14+OKv14_to_OIKv14)*OKv14;
a2 = C3Kv14_to_OKv14*C3Kv14 + OIKv14_to_OKv14*OIKv14;
dOKv14 = a2 - a1;

a1 = (CI0Kv14_to_C0Kv14+CI0Kv14_to_CI1Kv14)*CI0Kv14;
a2 = C0Kv14_to_CI0Kv14*C0Kv14 + CI1Kv14_to_CI0Kv14*CI1Kv14;
dCI0Kv14 = a2 - a1;

a1 = (CI1Kv14_to_CI2Kv14+CI1Kv14_to_C1Kv14+CI1Kv14_to_CI0Kv14)*CI1Kv14;
a2 = CI2Kv14_to_CI1Kv14*CI2Kv14 + C1Kv14_to_CI1Kv14*C1Kv14 + CI0Kv14_to_CI1Kv14*CI0Kv14;
dCI1Kv14 = a2 - a1;

a1 = (CI2Kv14_to_CI3Kv14+CI2Kv14_to_C2Kv14+CI2Kv14_to_CI1Kv14)*CI2Kv14;
a2 = CI3Kv14_to_CI2Kv14*CI3Kv14 + C2Kv14_to_CI2Kv14*C2Kv14 + CI1Kv14_to_CI2Kv14*CI1Kv14;
dCI2Kv14 = a2 - a1;

a1 = (CI3Kv14_to_OIKv14+CI3Kv14_to_C3Kv14+CI3Kv14_to_CI2Kv14)*CI3Kv14;
a2 = OIKv14_to_CI3Kv14*OIKv14 + C3Kv14_to_CI3Kv14*C3Kv14 + CI2Kv14_to_CI3Kv14*CI2Kv14;
dCI3Kv14 = a2 - a1;

a1 = (OIKv14_to_OKv14+OIKv14_to_CI3Kv14)*OIKv14;
a2 = OKv14_to_OIKv14*OKv14 + CI3Kv14_to_OIKv14*CI3Kv14;
dOIKv14 = a2 - a1;

                                    %   Assign global state derviative to derivative array
if Vclamp_flag
    dy(Nstates_CaRU+1) = 0;
else    
    dy(Nstates_CaRU+1) = dV;
end
dy(Nstates_CaRU+2) = dP_UO;%dmNa;
dy(Nstates_CaRU+3) = dP_UIF;%dhNa;
dy(Nstates_CaRU+4) = dP_UC2;%djNa;
dy(Nstates_CaRU+5) = dNai;
dy(Nstates_CaRU+6) = dKi;
dy(Nstates_CaRU+7) = dCai;
dy(Nstates_CaRU+8) = dCaSR;
dy(Nstates_CaRU+9) = dLTRPNCa;
dy(Nstates_CaRU+10) = dHTRPNCa;
dy(Nstates_CaRU+11) = dxKs;
dy(Nstates_CaRU+12) = dC0Kv43;
dy(Nstates_CaRU+13) = dC1Kv43;
dy(Nstates_CaRU+14) = dC2Kv43;
dy(Nstates_CaRU+15) = dC3Kv43;
dy(Nstates_CaRU+16) = dOKv43;
dy(Nstates_CaRU+17) = dCI0Kv43;
dy(Nstates_CaRU+18) = dCI1Kv43;
dy(Nstates_CaRU+19) = dCI2Kv43;
dy(Nstates_CaRU+20) = dCI3Kv43;
dy(Nstates_CaRU+21) = dOIKv43;
dy(Nstates_CaRU+22) = dC0Kv14;
dy(Nstates_CaRU+23) = dC1Kv14;
dy(Nstates_CaRU+24) = dC2Kv14;
dy(Nstates_CaRU+25) = dC3Kv14;
dy(Nstates_CaRU+26) = dOKv14;
dy(Nstates_CaRU+27) = dCI0Kv14;
dy(Nstates_CaRU+28) = dCI1Kv14;
dy(Nstates_CaRU+29) = dCI2Kv14;
dy(Nstates_CaRU+30) = dCI3Kv14;
dy(Nstates_CaRU+31) = dOIKv14;
dy(Nstates_CaRU+32) = dC1Herg;
dy(Nstates_CaRU+33) = dC2Herg;
dy(Nstates_CaRU+34) = dC3Herg;
dy(Nstates_CaRU+35) = dOHerg;
dy(Nstates_CaRU+36) = dIHerg;
dy(Nstates_CaRU+37) = dP_UC3;%dhNal;
dy(Nstates_CaRU+38) = dOIKv43Ca;
dy(Nstates_CaRU+39) = dP_UIM1;
dy(Nstates_CaRU+40) = dP_UC1;
dy(Nstates_CaRU+41) = dP_IC3;
dy(Nstates_CaRU+42) = dP_IC2;
dy(Nstates_CaRU+43) = dP_UIM2;
dy(Nstates_CaRU+44) = dP_LC3;
dy(Nstates_CaRU+45) = dP_LC2;
dy(Nstates_CaRU+46) = dP_LC1;
dy(Nstates_CaRU+47) = dP_LO;

if t > T_array(Counter)                %   Roughly eliminates data from rejected time steps
    Counter = Counter + 1;  
end
T_array(Counter) = t;                  %   Store currents/fluxes for output
ICaL_array(Counter) = ICaL;
JCaL_array(Counter) = JCaL;
JRyR_array(Counter) = Jup;
CaSSavg_array(Counter) = CaSSavg;
IKs_array(Counter,1) = (IKs);





Loading data, please wait...