% 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