- Author:
- Francesca Simone <francesca.simone01@universitadipavia.it>
- Date:
- 2026-01-14 15:17:24+01:00
- Desc:
- Add new CELLML model
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/d92/rawfile/83903a0ad0a326782c7f35c4d51e9cea5dc6e90a/Simone_Trancuccio2025.m
% Software implementation of the Simone_Trancuccio2025 model of the action potential
% of human induced pluripotent stem cell-derived cardiomyocytes.
%
% This software is provided for NON-COMMERCIAL USE ONLY
%
% The Simone_Trancuccio2025 model is described in the article
% " A Novel Computational Model of Human iPSC-Derived Ventricular Myocytes with
% Improved L-Type Calcium Current for Application to Timothy Syndrome"
% Francesca Simone, Alessandro Trancuccio, Jaroslaw Karol Sochacki, Celia
% Martinez Prieto, Silvia G. Priori, Luca F. Pavarino, and Demetrio J.
% Santiago.
%
% Please cite the above paper when using this model.
%
% Email: francesca.simone01@universitadipavia.it /alessandro.trancuccio@gmail.com
%
function [output] = Simone_Trancuccio2025(time, Y, stimFlag, odeflag,i_stim_frequency,fractional_TS)
%-------------------------------------------------------------------------------
% State variables
%-------------------------------------------------------------------------------
% 1: Vm (volt) (in Membrane)
% 2: Ca_SR (millimolar) (in calcium_dynamics)
% 3: Cai (millimolar) (in calcium_dynamics)
% NOT USED 4: g (dimensionless) (in calcium_dynamics)
% 5: Ok (dimensionless) (VDI)
% 6: I1k (dimensionless) (VDI)
% 7: I2k (dimensionless) (VDI)
% 8: Ck (dimensionless) (VDI)
% 9: I1Cak (dimensionless) (CDI)
% 10: I2Cak (dimensionless) (CDI)
% 11: CCak (dimensionless) (CDI)
% 12: nca (dimensionless) (n Gate: regola la transizione)
% 13: jnca (dimensionless) (n Gate: regola la transizione)
% 14: Xr1 (dimensionless) (in i_Kr_Xr1_gate)
% 15: Xr2 (dimensionless) (in i_Kr_Xr2_gate)
% 16: Xs (dimensionless) (in i_Ks_Xs_gate)
% 17: h (dimensionless) (in i_Na_h_gate)
% 18: j (dimensionless) (in i_Na_j_gate)
% 19: m (dimensionless) (in i_Na_m_gate)
% 20: Xf (dimensionless) (in i_f_Xf_gate)
% 21: q (dimensionless) (in i_to_q_gate)
% 22: r (dimensionless) (in i_to_r_gate)
% 23: Nai (millimolar) (in sodium_dynamics)
% 24: m_L (dimensionless) (in i_NaL_m_gate)
% 25: h_L (dimensionless) (in i_NaL_h_gate)
% 26: RyRa (dimensionless) (in calcium_dynamics)
% 27: RyRo (dimensionless) (in calcium_dynamics)
% 28: RyRc (dimensionless) (in calcium_dynamics)
% 29: Ok (dimensionless) (VDI)
% 30: I1k (dimensionless) (VDI)
% 31: I2k (dimensionless) (VDI)
% 31: Ck (dimensionless) (VDI)
% 33: I1Cak (dimensionless) (CDI)
% 34: I2Cak (dimensionless) (CDI)
% 35: CCak (dimensionless) (CDI)
% 36: nca (dimensionless) (n Gate: regola la transizione)
% 37: jnca (dimensionless) (n Gate: regola la transizione)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Constants
F = 96485.3415; % coulomb_per_mole (in model_parameters)
R = 8.314472; % joule_per_mole_kelvin (in model_parameters)
T = 310.0; % kelvin (in model_parameters) %37°C
vffrt = Y(1)*F*F/(R*T);
vfrt = Y(1)*F/(R*T);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Cell geometry
V_SR = 583.73; % micrometre_cube (in model_parameters)
Vc = 8800.0; % micrometre_cube (in model_parameters)
Cm = 9.87109e-11; % farad (in model_parameters)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Extracellular concentrations
Nao = 151.0; % millimolar (in model_parameters)
Ko = 5.4; % millimolar (in model_parameters)
Cao =1.8; % millimolar (in model_parameters)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Intracellular concentrations
% Naio = 10 mM Y(18)
Ki = 150.0; % millimolar (in model_parameters)
% Cai = 0.0002 mM Y(3)
% caSR = 0.3 mM Y(2)
% time (second)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Nernst potential
E_Na = R*T/F*log(Nao/Y(23));
E_Ca = 0.5*R*T/F*log(Cao/Y(3));
E_K = R*T/F*log(Ko/Ki);
PkNa = 0.03; % dimensionless (in electric_potentials)
E_Ks = R*T/F*log((Ko+PkNa*Nao)/(Ki+PkNa*Y(23)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INa adapted from DOI:10.3389/fphys.2018.00080
g_Na = 1.2176*6447.1896;
i_Na = g_Na*Y(19)^3.0*Y(17)*Y(18)*(Y(1) - E_Na);
m_inf = 1 / (1 + exp((Y(1)*1000 + 39)/-11.2));
tau_m = (0.00001 + 0.00013*exp(-((Y(1)*1000 + 48)/15)^2) + 0.000045 / (1 + exp((Y(1)*1000 + 42)/-5)));
dY(19, 1) = (m_inf-Y(19))/tau_m;
h_inf = 1 / (1 + exp((Y(1)*1000 + 66.5)/6.8));
tau_h = (0.00007 + 0.034 / (1 + exp((Y(1)*1000 + 41)/5.5) + exp(-(Y(1)*1000 + 41)/14)) + 0.0002 / (1 + exp(-(Y(1)*1000 + 79)/14)));
dY(17,1) = (h_inf-Y(17))/tau_h;
j_inf = h_inf;
tau_j = 10*(0.0007 + 0.15 / (1 + exp((Y(1)*1000 + 41)/5.5) + exp(-(Y(1)*1000 + 41)/14)) + 0.002 / (1 + exp(-(Y(1)*1000 + 79)/14)));
dY(18, 1) = (j_inf-Y(18))/tau_j;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INaL
myCoefTauM = 1;
tauINaL = 200; %ms
GNaLmax = 1.3021*2.3*7.5; %(S/F)
Vh_hLate = 87.61;
i_NaL = GNaLmax* Y(24)^(3)*Y(25)*(Y(1)-E_Na);
m_inf_L = 1/(1+exp(-(Y(1)*1000+42.85)/(5.264)));
alpha_m_L = 1/(1+exp((-60-Y(1)*1000)/5));
beta_m_L = 0.1/(1+exp((Y(1)*1000+35)/5))+0.1/(1+exp((Y(1)*1000-50)/200));
tau_m_L = 1/1000 * myCoefTauM*alpha_m_L*beta_m_L;
dY(24, 1) = (m_inf_L-Y(24))/tau_m_L;
h_inf_L = 1/(1+exp((Y(1)*1000+Vh_hLate)/(7.488)));
tau_h_L = 1/1000 * tauINaL;
dY(25, 1) = (h_inf_L-Y(25))/tau_h_L;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% If adapted from DOI:10.3389/fphys.2018.00080
g_f = 0.62395*22.2763088;
fNa = 0.37;
fK = 1 - fNa;
i_fK = fK*g_f*Y(20)*(Y(1) - E_K);
i_fNa = fNa*g_f*Y(20)*(Y(1) - E_Na);
i_f = i_fK + i_fNa;
Xf_infinity = 1.0/(1.0 + exp((Y(1)*1000 + 69)/8));
tau_Xf = 5600 / (1 + exp((Y(1)*1000 + 65)/7) + exp(-(Y(1)*1000 + 65)/19));
dY(20, 1) = 1000*(Xf_infinity-Y(20))/tau_Xf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ICaL
%%% ICaL parameters %%%
bGCaL =1;
kCDI = 4.4381;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% up/down rates
%%% CDI on/off
undo_CDI = 0;
% bn: ICaL CDI block (VDI-only ICaL)
% - bn=1* -> no CDI block
% - bn=0 -> CDI total block
bn=1;
r_down = (1e-1);
r_down_TS = (1e-1);
r_up = bn*(r_down*Y(12)/(1-Y(12)))*(1-undo_CDI);
r_up_TS=bn*(r_down*Y(36)/(1-Y(36)))*(1-undo_CDI);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% n gate -> used for compute nca
jncass = 1.0/(1.0+exp((Y(1)*1000+19.58+25)/3.696));
tjnca = 1.1286;
dY(13,1) =(jncass-Y(13))/tjnca;
%%n gate -> used for compute nca TS
jncass = 1.0/(1.0+exp((Y(1)*1000+19.58+25)/3.696));
tjnca = 1.1286;
dY(37,1) =(jncass-Y(37))/tjnca;
Kmn = 5.6e-05;
k2n = 1000;
km2n = 150*Y(13);
anca = (1-Y(12))/(1+Kmn/Y(3))^4.0;
dY(12,1)=bn*(anca*k2n-Y(12)*km2n);
%%TS
km2nTS = 150*Y(37);
ancaTS=(1-Y(36))/(1+Kmn/Y(3))^4.0;
dY(36,1)=bn*(ancaTS*k2n-Y(36)*km2nTS);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Activation (d) (C<->0)
dss = 1.0/(1.0+exp((-(Y(1)*1000+4.0918))/5.8312));
td = (0.000262+9.2e-5/( exp( -0.050025*(Y(1)*1000+6.3908))+exp(0.10331*(Y(1)*1000+14.5728))) );
alpha =dss/td;
beta =(1-dss) / td;
% Activation (d) TS
dss = 1.0/(1.0+exp((-(Y(1)*1000+4.0918))/5.8312));
td = (0.000262+9.2e-5/( exp( -0.050025*(Y(1)*1000+6.3908))+exp(0.10331*(Y(1)*1000+14.5728))) );
alpha_TS =dss/td;
beta_TS =(1-dss) / td;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recovery (jca)
jcass_new= 1.0/(1.0+exp((Y(1)*1000+19.58)/3.696));
jcass_VD = jcass_new;
jcass_CD = jcass_new;
tjca_new = 35 + 350*exp(-(Y(1)*1000-(-20))^2/(2*100));
tjca_VD = tjca_new;
tjca_CD = tjca_new;
% psi and omega rates
psi_VD=jcass_VD/tjca_VD;
psi_CD=jcass_CD/tjca_CD;
omega_VD=(1-jcass_VD)/tjca_VD;
omega_CD=(1-jcass_CD)/tjca_CD;
psi_VD_TS=jcass_VD/tjca_VD;
psi_CD_TS=jcass_CD/tjca_CD;
omega_VD_TS=0.1*(1-jcass_VD)/tjca_VD;
omega_CD_TS=0.1*(1-jcass_CD)/tjca_CD;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fact Inactivation (f1) (0 <-> I_1)
f1ss_0 = (0.092141 / (1.0+exp((Y(1)*1000+20.2666)/3.7342)) + 0.017445);
tf1_0 = 1*(0.092321+ 0.00094./ (0.005431*exp((Y(1)*1000+22.3472)/(-46.5268))+0.005082*exp((Y(1)*1000+32.8386)/11.3505)));
gamma_VD = (1-f1ss_0)/ tf1_0;
delta_VD = f1ss_0 / tf1_0;
gamma_CD=gamma_VD*kCDI;
delta_CD=delta_VD*kCDI;
%TS
gamma_VD_TS = 0.06*(1-f1ss_0)/ tf1_0;
delta_VD_TS = f1ss_0 / tf1_0;
gamma_CD_TS=gamma_VD*kCDI;
delta_CD_TS=delta_VD*kCDI;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slow Inactivation (f2)
tf2_new = 1*(106.9697 + 0./ (0.004587*exp((Y(1)*1000+4.957)/(-90.6055))+0.003639*exp((Y(1)*1000+4.6196)/4.0157)));
tf2_VD = tf2_new;
tf2_CD = tf2_VD/kCDI;
%TS
tf2_new_TS = 1*(106.9697 + 0./ (0.004587*exp((Y(1)*1000+4.957)/(-90.6055))+0.003639*exp((Y(1)*1000+4.6196)/4.0157)));
tf2_VD_TS = tf2_new_TS;
tf2_CD_TS = tf2_VD_TS/kCDI;
% Reversibility
theta_VD = alpha*gamma_VD*psi_VD/tf2_VD/(alpha*gamma_VD*psi_VD+beta*delta_VD*omega_VD);
theta_CD = alpha*gamma_CD*psi_CD/tf2_CD/(alpha*gamma_CD*psi_CD+beta*delta_CD*omega_CD);
eta_VD=1/tf2_VD - theta_VD;%%
eta_CD=1/tf2_CD - theta_CD;%%
%TS
theta_VD_TS = alpha_TS*gamma_VD_TS*psi_VD_TS/tf2_VD_TS/(alpha_TS*gamma_VD_TS*psi_VD_TS+beta_TS*delta_VD_TS*omega_VD_TS);
theta_CD_TS = alpha_TS*gamma_CD_TS*psi_CD_TS/tf2_CD_TS/(alpha_TS*gamma_CD_TS*psi_CD_TS+beta_TS*delta_CD_TS*omega_CD_TS);
eta_VD_TS=1/tf2_VD - theta_VD;%%
eta_CD_TS=1/tf2_CD - theta_CD;%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Driving Forces
PhiCaL=4.0*vffrt*(1.2*Y(3)*exp(2.0*vfrt)-0.341*Cao)/(exp(2*vfrt)-1.0);
PhiCaNa=1.0*vffrt*(0.75*Y(23)*exp(1.0*vfrt)-0.75*Nao)/(exp(1.0*vfrt)-1.0);
PhiCaK=1.0*vffrt*(0.75*Ki*exp(1.0*vfrt)-0.75*Ko)/(exp(1.0*vfrt)-1.0);
PCa=0.87069*0.000147;
PCaNa=0.00125*PCa;
PCaK=3.574e-4*PCa;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Markov Model: VDI states
OCaK = 1-Y(11)-Y(9)-Y(10)-Y(8)-Y(6)-Y(7)-Y(5);
dY(5,1) = alpha*Y(8) + delta_VD*Y(6) - (beta+gamma_VD)*Y(5) - r_up*Y(5) + r_down*OCaK;
dY(7,1) = eta_VD*Y(6) + omega_VD*Y(8) - (theta_VD+psi_VD)*Y(7) - r_up*Y(7) + r_down*Y(10);
dY(6,1) = theta_VD*Y(7) + gamma_VD*Y(5) - (eta_VD+delta_VD)*Y(6) - r_up*Y(6) + r_down*Y(9);
dY(8,1) = beta*Y(5) + psi_VD*Y(7) - (omega_VD+alpha)*Y(8) - r_up*Y(8) + r_down*Y(11);
%TS
OCaK_TS = 1-Y(35)-Y(33)-Y(34)-Y(32)-Y(30)-Y(31)-Y(29);
dY(29,1) = alpha_TS*Y(32) + delta_VD_TS*Y(30) - (beta_TS+gamma_VD_TS)*Y(29) - r_up_TS*Y(29) + r_down_TS*OCaK_TS;
dY(31,1) = eta_VD_TS*Y(30) + omega_VD_TS*Y(32) - (theta_VD_TS+psi_VD_TS)*Y(31) - r_up_TS*Y(31) + r_down_TS*Y(34);
dY(30,1) = theta_VD_TS*Y(31) + gamma_VD_TS*Y(29) - (eta_VD_TS+delta_VD_TS)*Y(30) - r_up_TS*Y(30) + r_down_TS*Y(33);
dY(32,1) = beta_TS*Y(29) + psi_VD_TS*Y(31) - (omega_VD_TS+alpha_TS)*Y(32) - r_up_TS*Y(32) + r_down_TS*Y(35);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Markov Model: CDI states
dY(10,1) = eta_CD*Y(9) + omega_CD*Y(11) - (theta_CD+psi_CD)*Y(10)+ r_up*Y(7) - r_down*Y(10);
dY(9,1) = theta_CD*Y(10) + gamma_CD*OCaK - (eta_CD+delta_CD)*Y(9) + r_up*Y(6) - r_down*Y(9);
dY(11,1) = beta*OCaK + psi_CD*Y(10) - (omega_CD+alpha)*Y(11) + r_up*Y(8) - r_down*Y(11);
%TS
dY(34,1) = eta_CD_TS*Y(33) + omega_CD_TS*Y(35) - (theta_CD_TS+psi_CD_TS)*Y(34)+ r_up_TS*Y(31) - r_down_TS*Y(34);
dY(33,1) = theta_CD_TS*Y(34) + gamma_CD_TS*OCaK_TS - (eta_CD_TS+delta_CD_TS)*Y(33) + r_up_TS*Y(30) - r_down_TS*Y(33);
dY(35,1) = beta_TS*OCaK_TS + psi_CD_TS*Y(34) - (omega_CD_TS+alpha_TS)*Y(35) + r_up_TS*Y(32) - r_down_TS*Y(35);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reversibility
check_reversibility=0;
if check_reversibility > 0
revTol = check_reversibility;
rev1=abs(alpha*gamma_VD*eta_VD*psi_VD-beta*delta_VD*theta_VD*omega_VD);
rev2=abs(alpha*gamma_VDp*eta_VDp*psi_VDp-beta*delta_VDp*theta_VDp*omega_VDp);
rev3=abs(alpha*gamma_CD*eta_CD*psi_CD-beta*delta_CD*theta_CD*omega_CD);
rev4=abs(alpha*gamma_CDp*eta_CDp*psi_CDp-beta*delta_CDp*theta_CDp*omega_CDp);
if rev1>revTol || rev2>revTol || rev3>revTol || rev4>revTol
disp('REVERSIBILITY FAILED');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICaL ICaNa ICaK currents
Ib=1;
ICaL_VD = Ib*PCa *PhiCaL *Y(5);
ICaL_CD = Ib*PCa *PhiCaL *OCaK;
ICaNa_VD = Ib*PCaNa *PhiCaNa *Y(5);
ICaNa_CD = Ib*PCaNa *PhiCaNa *OCaK;
ICaK_VD = Ib*PCaK *PhiCaK *Y(5);
ICaK_CD = Ib*PCaK *PhiCaK *OCaK;
%TS
ICaL_VD_TS = Ib*PCa *PhiCaL *Y(29);
ICaL_CD_TS = Ib*PCa *PhiCaL *OCaK_TS;
ICaNa_VD_TS = Ib*PCaNa *PhiCaNa *Y(29);
ICaNa_CD_TS = Ib*PCaNa *PhiCaNa *OCaK_TS;
ICaK_VD_TS = Ib*PCaK *PhiCaK *Y(29);
ICaK_CD_TS = Ib*PCaK *PhiCaK *OCaK_TS;
% ICaL VD vs CD & ICaL p vs np +TS
ICaLnp = ICaL_VD + ICaL_CD;
ICaLnp_TS = ICaL_VD_TS + ICaL_CD_TS;
ICaNanp = ICaNa_VD + ICaNa_CD;
ICaNanp_TS = ICaNa_VD_TS + ICaNa_CD_TS;
ICaKnp = ICaK_VD + ICaK_CD;
ICaKnp_TS = ICaK_VD_TS + ICaK_CD_TS;
ICaL = (ICaLnp)*bGCaL;
ICaL_TS = (ICaLnp_TS)*bGCaL;
ICaNa = ( ICaNanp)*bGCaL;
ICaNa_TS = ( ICaNanp_TS)*bGCaL;
ICaK = (ICaKnp)*bGCaL;
ICaK_TS = (ICaKnp_TS)*bGCaL;
i_CaL=(1-fractional_TS)*(ICaL+ICaNa+ICaK)+(fractional_TS)*(ICaL_TS+ICaNa_TS+ICaK_TS);
% ICaL conductance
gICaL = ICaL/PhiCaL;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%% Ito
g_to = 1.2913*29.9038; % S_per_F (in i_to)
i_to = g_to*(Y(1)-E_K)*Y(21)*Y(22);
q_inf = 1.0/(1.0+exp((Y(1)*1000.0+53.0)/13.0));
tau_q = (6.06+39.102/(0.57*exp(-0.08*(Y(1)*1000.0+44.0))+0.065*exp(0.1*(Y(1)*1000.0+45.93))))/1000.0;
dY(21, 1) = (q_inf-Y(21))/tau_q;
r_inf = 1.0/(1.0+exp(-(Y(1)*1000.0-22.3)/18.75));
tau_r = (2.75352+14.40516/(1.037*exp(0.09*(Y(1)*1000.0+30.61))+0.369*exp(-0.12*(Y(1)*1000.0+23.84))))/1000.0;
dY(22, 1) = (r_inf-Y(22))/tau_r;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IKs
g_Ks = 1.0726*2.041; % S_per_F (in i_Ks)
i_Ks =g_Ks*(Y(1)-E_Ks)*Y(16)^2.0*(1.0+0.6/(1.0+(3.8*0.00001/Y(3))^1.4));
Xs_infinity = 1.0/(1.0+exp((-Y(1)*1000.0-20.0)/16.0));
alpha_Xs = 1100.0/sqrt(1.0+exp((-10.0-Y(1)*1000.0)/6.0));
beta_Xs = 1.0/(1.0+exp((-60.0+Y(1)*1000.0)/20.0));
tau_Xs = 1.0*alpha_Xs*beta_Xs/1000.0;
dY(16, 1) = (Xs_infinity-Y(16))/tau_Xs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IKr
L0 = 0.025; % dimensionless (in i_Kr_Xr1_gate)
Q = 2.3; % dimensionless (in i_Kr_Xr1_gate)
g_Kr = 1.9913*29.8667; % S_per_F (in i_Kr)
i_Kr = g_Kr*(Y(1)-E_K)*Y(14)*Y(15)*sqrt(Ko/5.4);
V_half = 1000.0*(-R*T/(F*Q)*log((1.0+Cao/2.6)^4.0/(L0*(1.0+Cao/0.58)^4.0))-0.019);
Xr1_inf = 1.0/(1.0+exp((V_half-Y(1)*1000.0)/4.9));
alpha_Xr1 = 450.0/(1.0+exp((-45.0-Y(1)*1000.0)/10.0));
beta_Xr1 = 6.0/(1.0+exp((30.0+Y(1)*1000.0)/11.5));
tau_Xr1 = 1.0*alpha_Xr1*beta_Xr1/1000.0;
dY(14, 1) = (Xr1_inf-Y(14))/tau_Xr1;
Xr2_infinity = 1.0/(1.0+exp((Y(1)*1000.0+88.0)/50.0));
alpha_Xr2 = 3.0/(1.0+exp((-60.0-Y(1)*1000.0)/20.0));
beta_Xr2 = 1.12/(1.0+exp((-60.0+Y(1)*1000.0)/20.0));
tau_Xr2 = 1.0*alpha_Xr2*beta_Xr2/1000.0;
dY(15, 1) = (Xr2_infinity-Y(15))/tau_Xr2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IK1
alpha_K1 = 3.91/(1.0+exp(0.5942*(Y(1)*1000.0-E_K*1000.0-200.0)));
beta_K1 = (-1.509*exp(0.0002*(Y(1)*1000.0-E_K*1000.0+100.0))+exp(0.5886*(Y(1)*1000.0-E_K*1000.0-10.0)))/(1.0+exp(0.4547*(Y(1)*1000.0-E_K*1000.0)));
XK1_inf = alpha_K1/(alpha_K1+beta_K1);
g_K1 = 0.50122*28.1492; % S_per_F (in i_K1)
i_K1 = g_K1*XK1_inf*(Y(1)-E_K)*sqrt(Ko/5.4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INaCa
KmCa = 1.38; % millimolar (in i_NaCa)
KmNai = 87.5; % millimolar (in i_NaCa)
Ksat = 0.1; % dimensionless (in i_NaCa)
gamma = 0.35; % dimensionless (in i_NaCa)
alpha1 = 2.16659;
kNaCa = 0.66821*6514.47574; % A_per_F (in i_NaCa)
i_NaCa = kNaCa*(exp(gamma*Y(1)*F/(R*T))*Y(23)^3.0*Cao-exp((gamma-1.0)*Y(1)*F/(R*T))*Nao^3.0*Y(3)*alpha1)/((KmNai^3.0+Nao^3.0)*(KmCa+Cao)*(1.0+Ksat*exp((gamma-1.0)*Y(1)*F/(R*T))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%% INaK
Km_K = 1.0; % millimolar (in i_NaK)
Km_Na = 40.0; % millimolar (in i_NaK)
PNaK = 1.3055*2.74240;% A_per_F (in i_NaK)
i_NaK = PNaK*Ko/(Ko+Km_K)*Y(23)/(Y(23)+Km_Na)/(1.0+0.1245*exp(-0.1*Y(1)*F/(R*T))+0.0353*exp(-Y(1)*F/(R*T)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% IpCa
KPCa = 0.0005; % millimolar (in i_PCa)
g_PCa = 1.5351*0.4125; % A_per_F (in i_PCa)
i_PCa = g_PCa*Y(3)/(Y(3)+KPCa);
%% Background currents
g_b_Na = 1.14; % S_per_F (in i_b_Na)
i_b_Na = g_b_Na*(Y(1)-E_Na);
g_b_Ca = 0.8727264; % S_per_F (in i_b_Ca)
i_b_Ca = g_b_Ca*(Y(1)-E_Ca);
%% Sarcoplasmic reticulum
VmaxUp = 1.4851*0.82205;
Kup = 0.52036* 4.40435e-4;
i_up = VmaxUp/(1.0+Kup^2.0/Y(3)^2.0);
V_leak = 4.48209e-4;
i_leak = (Y(2)-Y(3))*V_leak;
dY(4, 1) = 0;
% RyR
g_irel_max = 0.71915*55.808061;
RyRa1 = 0.05169;
RyRa2 = 0.050001;
RyRahalf = 0.02632;
RyRohalf = 0.00944;
RyRchalf = 0.00167;
RyRSRCass = (1 - 1/(1 + exp((Y(2)-0.3)/0.1)));
i_rel = g_irel_max*RyRSRCass*Y(27)*Y(28)*(Y(2)-Y(3));
RyRainfss = RyRa1-RyRa2/(1 + exp((1000*Y(3)-(RyRahalf))/0.0082));
RyRtauadapt = 1; %s
dY(26,1) = (RyRainfss- Y(26))/RyRtauadapt;
RyRoinfss = (1 - 1/(1 + exp((1000*Y(3)-(Y(26)+ RyRohalf))/0.003)));
if (RyRoinfss>= Y(27))
RyRtauact = 18.75e-3; %s
else
RyRtauact = 0.1*18.75e-3; %s
end
dY(27,1) = (RyRoinfss- Y(27))/(RyRtauact);
RyRcinfss = (1/(1 + exp((1000*Y(3)-(Y(26)+RyRchalf))/0.001)));
if (RyRcinfss>= Y(28))
RyRtauinact = 2*87.5e-3; %s
else
RyRtauinact = 87.5e-3; %s
end
dY(28,1) = (RyRcinfss- Y(28))/(RyRtauinact);
%% Ca2+ buffering
Buf_C = 0.25; % millimolar (in calcium_dynamics)
Buf_SR = 10.0; % millimolar (in calcium_dynamics)
Kbuf_C = 0.001; % millimolar (in calcium_dynamics)
Kbuf_SR = 0.3; % millimolar (in calcium_dynamics)
Cai_bufc = 1.0/(1.0+Buf_C*Kbuf_C/(Y(3)+Kbuf_C)^2.0);
Ca_SR_bufSR = 1.0/(1.0+Buf_SR*Kbuf_SR/(Y(2)+Kbuf_SR)^2.0);
%% Ionic concentrations
%Nai
dY(23, 1) = -Cm*(i_Na+i_NaL+i_b_Na+3.0*i_NaK+3.0*i_NaCa+i_fNa)/(F*Vc*1.0e-18);
%Cai
dY(3, 1) = Cai_bufc*(i_leak-i_up+i_rel-(i_CaL+i_b_Ca+i_PCa-2.0*i_NaCa)*Cm/(2.0*Vc*F*1.0e-18));
%caSR
dY(2, 1) = Ca_SR_bufSR*Vc/V_SR*(i_up-(i_rel+i_leak));
%% Stimulation
i_stim_Amplitude = 5.5e-10;%7.5e-10; % ampere (in stim_mode)
i_stim_End = 1000.0; % second (in stim_mode)
i_stim_PulseDuration = 0.005; % second (in stim_mode)
i_stim_Start = 0; % second (in stim_mode) %0.1
stim_flag = stimFlag; % dimensionless (in stim_mode)
i_stim_Period = 1/i_stim_frequency; % second
if stim_flag~=0 && stim_flag~=1
error('Simone_Trancuccio2025: wrong pacing! stimFlag can be only 0 (spontaneous) or 1 (paced)');
end
if ((time >= i_stim_Start) && (time <= i_stim_End) && (time-i_stim_Start-floor((time-i_stim_Start)/i_stim_Period)*i_stim_Period <= i_stim_PulseDuration))
i_stim = stim_flag*i_stim_Amplitude/Cm;
else
i_stim = 0.0;
end
%% Membrane potential
dY(1, 1) = -(i_K1+i_to+i_Kr+i_Ks+i_CaL+i_NaK+i_Na+i_NaL+i_NaCa+i_PCa+i_f+i_b_Na+i_b_Ca-i_stim);
%% Output variables
IK1 = i_K1;
Ito = i_to;
IKr = i_Kr;
IKs = i_Ks;
ICaL = i_CaL;
INaK = i_NaK;
INa = i_Na;
INaCa = i_NaCa;
IpCa = i_PCa;
If = i_f;
IbNa = i_b_Na;
IbCa = i_b_Ca;
Irel = i_rel;
Iup = i_up;
Ileak = i_leak;
Istim = i_stim;
INaL = i_NaL;
if odeflag ==1
output=dY;
else
output=[INa, If, ICaL, Ito, IKs, IKr, IK1, INaCa, INaK, IpCa, IbNa, IbCa, Irel, Istim, E_K, E_Na, INaL, Iup];
end
