clc, clear all, close all
%%%%%%%%%%%%%% Main File to run the Trovato2020 model %%%%%%%%%%%%%%%%%%%%%
% The following code simulates 1 beat at 1 Hz with                        %
% the human cardiac Purkinje Trovato2020 model                            %
% -                                                                       %  
% Model description and validation in:                                    %
% C. Trovato, E. Passini, N. Nagy, et al., Human Purkinje in              %
% silico model enables mechanistic investigations into automaticity       %
% and pro-arrhythmic abnormalities, Journal of Molecular and Cellular     %
% Cardiology (2019), https://doi.org/10.1016/j.yjmcc.2020.04.001 .        %
% [open-access]                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model=@modTrovato2020;
tsim= 1000; %Simulation time (ms)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initial condition from Steady State [1000 beats at 1Hz]
v= -86.550102957989600;
nai= 8.231857731510893;     nasl= 8.231553454361393;     nass= 8.231561539013950;
ki= 1.437673447706863e+02;  ksl= 1.437677771226899e+02;  kss= 1.437677533510394e+02;
cai= 4.360080908582633e-05; casl= 1.020101597047671e-04; cass= 1.017658486729359e-04;
cansr= 1.263525645841406;   cajsr= 1.248146625349512;    cacsr = 1.265185617721750;
m= 0.006341207769833;       hf= 0.788541761218318;       hs= 0.788474863764949;
j= 0.790412100577539;       hsp= 0.579594138981772;      jp= 0.790885181877794;
mL= 0.0;                    hL= 0.463478975643765;       hLp= 0.240123148015689;
a= 0.0;                     i= 0.649386868151536;        i2= 0.989963717273401;
d= 0.0;                     ff= 0.999999963501871;       fs= 0.926598485399264;
fcaf= 0.999999963493016;    fcas= 0.999834100785052;     jca= 0.999979200703676;
nca= 0.005470852996192;     ffp= 0.999999963494571;      fcafp= 0.999999963612862;
b= 0.0;                     g= 0.994211562437775;        xrf= 0.0;
xrs= 0.568856244015729;     xs1= 0.191294664752654;      xs2= 0.0;
xk1= 0.997085416662044;     y= 0.233014639857230;        CaMKt= 0.005060490773142;
u= 0.466232550741101;       Jrel1= 0;                    Jrel2= 0;

X0=[v      nai    nasl   nass   ki... 
    ksl    kss    cai    casl   cass...
    cansr  cajsr  cacsr  m      hf...
    hs     j      hsp    jp     mL...
    hL     hLp    a      i      i2...
    d      ff     fs     fcaf   fcas...
    jca    nca    ffp    fcafp  b...
    g      xrf    xrs    xs1    xs2...
    y      xk1    CaMKt  u      Jrel1...
    Jrel2]'; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Simulation and AP Plot
[t,y] = ode15s(model,[0 tsim],X0);
figure(); plot(t,y(:,1)); xlabel('Time (ms)'); ylabel('V_{m} (mV)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%