% 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
