driver.m

%% INITIALIZATION

clear variables;
close all;
clc;

%% SETTINGS

% for example:
settings.BCL = 800;
settings.TSim = 900000; % 15 min
%settings.TSim = 900000+300000; % 20 min

settings.NumStim = floor(settings.TSim / settings.BCL);  % Number of action potentials to calculate
settings.storeLast = settings.NumStim; % I only want the very last one of each series of stimulation

settings.pH_status = 1; % pH decrescente
%settings.pH_status = 0; % pH fisso

settings.saving_status = 0; % salva ogni minuto
%settings.saving_status = 0; % salva tutto

[StateControl,TiControl, Ns, pH, extra_var]= maleckar_main(settings);

names_sv = ["Ca_{rel}", "Ca_{up}", "F1", "F2", "O_{Calse}", "r", "s", "d_L", "f_{L1}", "f_{L2}", "Ca_c", "K_c", "Na_c", "n", "pa", "O", "O_C", "O_{TC}", "O_{TMgC}", "O_{TMgMg}", "Ca_d", "Ca_i", "K_i", "Na_i", "V", "h1", "h2", "m", "a_{ur}", "i_{ur}"];
units_sv = ["mM", "mM", "–", "–", "–", "–", "–", "–", "–", "–", "mM", "mM", "mM", "–", "–", "–", "–", "–", "–", "–", "mM", "mM", "mM", "mM", "mV", "–", "–", "–", "–", "–"];

%clc;

% Time vector
time = TiControl/1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

[eigenvalues, J] = calculate_eigenvalues_jacobian(settings);

cond(J)

figure('color','w');
spy(J,30,'b')
set(gca, 'FontSize', 17, 'LineWidth', 2);
title(sprintf('cond(J) = %.0d',cond(J)),'FontSize',25);



% Estrazione della parte reale degli autovalori
real_parts = real(eigenvalues);

% Calcolo del rapporto richiesto
max_real = max(abs(real_parts));
min_real = min(abs(real_parts));
ratio = max_real / min_real; % stiffness ratio




%% Imposed pH over time

figure;
if(settings.pH_status == 1)
    plot(Ns*settings.BCL/60000, pH,'Marker','*','Color',[1, 0.6, 0.6],'MarkerSize',2);
    grid on;
    title(sprintf('pH Decrease in phase 1A of myocardial ischemia'));
else
    plot(Ns*settings.BCL/60000, pH,'Marker','*','Color',[0.2, 1, 0.4],'MarkerSize',2);
    grid on;
    title(sprintf('pH over time'));
end
xlabel('Time [min]');
ylabel('pH [–]');
set(gca, 'FontSize', 17, 'LineWidth', 2);

%% PLOT WITH ALL VARIABLES

% Number of variables
num_vars = size(StateControl, 2); % variabili di stato di default

% Layout for subplots
num_rows = ceil(sqrt(num_vars));
num_cols = ceil(num_vars / num_rows);

% Define the number of distinct colors needed
num_colors = max(30, num_vars); % Ensure at least 30 colors

% Generate the colormap
color_map = hsv(num_colors);

% Define the figure size
figure_size = [100, 100, 1200, 800]; % [left, bottom, width, height]

% Create the figure with the specified size
figure('Position', figure_size);
% Loop through each variable and plot it
for i = 1:num_vars
    subplot(num_rows, num_cols, i);
    plot(time/60, StateControl(:, i), 'LineWidth', 2, 'Color', color_map(mod(i-1, num_colors) + 1, :));
    grid on;
    title(sprintf('SV%d: %s', i, names_sv{i}));
    xlabel('Time [min]');
    ylabel(sprintf('%s [%s]', names_sv{i}, units_sv{i}));
    set(gca, 'FontSize', 12, 'LineWidth', 1.5);
end 

clear i;

if (settings.pH_status == 1)
    sgtitle(sprintf('State Variables @ pH DECREASING from %.2f to %.2f in %d min',pH(1),pH(end),time(end)/60));
elseif (settings.pH_status == 0)
        sgtitle(sprintf('State Variables @ pH = %.2f FIXED over %d min',pH(1),time(end)/60));
end

%% IONIC CONCENTRATIONS 

figure;
subplot(2,1,1)
data1 = StateControl(:, 24);
data2 = StateControl(:, 23);
x = time/60;
yyaxis left
plot(x, data1, 'LineWidth', 2);
xticks(0:1:15)
ylabel(sprintf('%s [%s]', names_sv(24), units_sv(24)));
hold on
grid on
yyaxis right
plot(x, data2, 'LineWidth', 2);
ylabel(sprintf('%s [%s]', names_sv(23), units_sv(23)));
xlabel('Time [min]'); % Aggiungi una label all'asse x
set(gca, 'FontSize',12,'LineWidth',1.5)
subplot(2,1,2)
chosen = 22;
plot(x, StateControl(:, chosen),'k','LineWidth',2);
xticks(0:1:15)
grid on
xlabel('Time [min]');
ylabel(sprintf('%s [%s]', names_sv(chosen), units_sv(chosen)));
set(gca, 'FontSize',12,'LineWidth',1.5)
sgtitle(sprintf('Relevant ionic species: intracellular concentration trends'));

clear x data1 data2;

%% SINGLE PLOT WITH USER-CHOSEN VARIABLE

disp('This is the list of State Variables:');
disp('1: Ca_rel (millimolar) (in Ca_handling_by_the_SR)');
disp('2: Ca_up (millimolar) (in Ca_handling_by_the_SR)');
disp('3: F1 (dimensionless) (in Ca_handling_by_the_SR)');
disp('4: F2 (dimensionless) (in Ca_handling_by_the_SR)');
disp('5: O_Calse (dimensionless) (in Ca_handling_by_the_SR)');
disp('6: r (dimensionless) (in Ca_independent_transient_outward_K_current_r_gate)');
disp('7: s (dimensionless) (in Ca_independent_transient_outward_K_current_s_gate)');
disp('8: d_L (dimensionless) (in L_type_Ca_channel_d_L_gate)');
disp('9: f_L1 (dimensionless) (in L_type_Ca_channel_f_L1_gate)');
disp('10: f_L2 (dimensionless) (in L_type_Ca_channel_f_L2_gate)');
disp('11: Ca_c (millimolar) (in cleft_space_ion_concentrations)');
disp('12: K_c (millimolar) (in cleft_space_ion_concentrations)');
disp('13: Na_c (millimolar) (in cleft_space_ion_concentrations)');
disp('14: n (dimensionless) (in delayed_rectifier_K_currents_n_gate)');
disp('15: pa (dimensionless) (in delayed_rectifier_K_currents_pa_gate)');
disp('16: O (dimensionless) (in intracellular_Ca_buffering)');
disp('17: O_C (dimensionless) (in intracellular_Ca_buffering)');
disp('18: O_TC (dimensionless) (in intracellular_Ca_buffering)');
disp('19: O_TMgC (dimensionless) (in intracellular_Ca_buffering)');
disp('20: O_TMgMg (dimensionless) (in intracellular_Ca_buffering)');
disp('21: Ca_d (millimolar) (in intracellular_ion_concentrations)');
disp('22: Ca_i (millimolar) (in intracellular_ion_concentrations)');
disp('23: K_i (millimolar) (in intracellular_ion_concentrations)');
disp('24: Na_i (millimolar) (in intracellular_ion_concentrations)');
disp('25: V (millivolt) (in membrane)');
disp('26: h1 (dimensionless) (in sodium_current_h1_gate)');
disp('27: h2 (dimensionless) (in sodium_current_h2_gate)');
disp('28: m (dimensionless) (in sodium_current_m_gate)');
disp('29: a_ur (dimensionless) (in ultra_rapid_K_current_aur_gate)');
disp('30: i_ur (dimensionless) (in ultra_rapid_K_current_iur_gate)');
%chosen = input('Which variable do you want to plot? Enter the corresponding number:\n');
chosen = 25; % fixed on V

while ~isnumeric(chosen) || chosen < 1 || chosen > 30
    disp('Error: You must insert a valid number (between 1 and 30).');
    chosen = input('Which variable do you want to plot? Enter the corresponding number:\n');
end

clc;

figure;
plot(time/60, StateControl(:, chosen),'k','LineWidth',2);
grid on
title(sprintf('SV%d: %s', chosen, names_sv(chosen)));
xlabel('Time [min]');
ylabel(sprintf('%s [%s]', names_sv(chosen), units_sv(chosen)));
set(gca, 'FontSize',12,'LineWidth',1.5)

%% EXTRA (CONTROL) VARIABLES (CURRENTS & CONDUCTANCES)

%% currents
i_Na = extra_var(:,1);
i_Ca_L = extra_var(:,2);
i_NaCa = extra_var(:,3);
...


i_imp = [i_Na, i_Ca_L, i_NaCa];
i_tit = {'i_{Na}', 'i_{CaL}', 'i_{NaCa}', 'i_{NaK}','i_{K_{s}}','i_{K_{r}'};

figure_size = [1,629,1710,361]; % [left, bottom, width, height]

% Create the figure with the specified size
figure('Position',figure_size,'Color','w');
for el = 1:size(i_imp,2)
    subplot(1,3,el)
    if el==1
    color =[0.2,0.6,0.8];
    elseif el==2
    color =[1,0.8,0];
    elseif el==3
    color =[0,0.6,0.2];
    end
    plot(time/60, i_imp(:,el),'Color',color,'LineWidth',2.5);
    xticks(0:1:15)
    grid on
    xlabel('Time [min]');
    ylabel('Current Density [pA/pF]');
    title(i_tit(el));
    set(gca, 'FontSize',17,'LineWidth',1.8)
end

%%

g_imp = [i_Na./(StateControl(:,25)-extra_var(:,9)), i_Ca_L./(StateControl(:,25)-extra_var(:,11))];
g_tit = {'g_{Na}', 'g_{CaL}', 'g_{NaCa}', 'g_{NaK}','g_{K_{s}}','g_{K_{r}'};

figure_size = [1,629,1710,361]; % [left, bottom, width, height]

% Create the figure with the specified size
figure('Position', figure_size,'Color','w');
for el = 1:size(g_imp,2)
    subplot(1,3,el)
    if el==1
    color =[0.2,0.6,0.8];
    elseif el==2
    color =[1,0.8,0];
    elseif el==3
    color =[0,0.6,0.2];
    end
    plot(time/60, g_imp(:,el),'Color',color,'LineWidth',2.5);
    xticks(0:1:15)
    grid on
    xlabel('Time [min]');
    ylabel('Conductance Density [pA/pF/mV]');
    title(g_tit(el));
    set(gca, 'FontSize',17,'LineWidth',1.8)
end

%% (conductances) --> not very significant ...

g_Na_max = extra_var(:,7);
g_Ca_L_max = extra_var(:,8);

figure;
plot(time/60, g_Na_max,'k','LineWidth',2);
grid on
xlabel('Time [min]');
ylabel('Current [nL/s]');
title('g_{Na_{max}}');
set(gca, 'FontSize',12,'LineWidth',1.5)

figure;
plot(time/60, g_Ca_L_max,'k','LineWidth',2);
grid on
xlabel('Time [min]');
ylabel('Current [nS]');
title('g_{CaL_{max}}');
set(gca, 'FontSize',12,'LineWidth',1.5)


%% PA a confronto 

chosen = 25;
V_div = [];
time_div = [];

BCL = settings.BCL / 1000; %[s]
Nstim = 60000 / settings.BCL; % [–]

% In queste variabili vengono salvati i potenziali di interesse e i relativi tempi
V_div_cell = cell(1, settings.TSim / 60000 +1);  % Inizializza il cell array con dimensioni corrispondenti al numero di iterazioni
time_div_cell = cell(1, settings.TSim / 60000 +1);  % Inizializza il cell array con dimensioni corrispondenti al numero di iterazioni

% Initialize the video writer
videoFileName = 'plot_video.mp4';
v = VideoWriter(videoFileName,  'MPEG-4');
v.FrameRate = 2; % Adjust frame rate as needed
v.Quality = 100;
open(v);

figure('Color', 'w'); % HD resolution
num_colors = max(1, (settings.TSim / 60000 +1));
color_map = hsv(num_colors);

for ord = 1:(settings.TSim / 60000 +1)
    if (ord == 1)
        for indice = 1:length(time)

            if (0<= time(indice) && time(indice) < BCL)
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    elseif (ord == 2)
        for indice = 1:length(time)
            if ((BCL < time(indice) && time(indice) <= BCL*Nstim))
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    elseif (ord > 2)
        for indice = 1:length(time)
            if ((ord-2)*BCL*Nstim < time(indice) && time(indice) <= (ord-1)*BCL*Nstim)
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    end


    % Aggiungi V_div al cell array
    V_div_cell{ord} = V_div;
    time_div_cell{ord} = time_div;

    % Plot una volta per ogni ord
    plot(1000*(time_div-time_div(1)),V_div, 'LineWidth', 2, 'Color', color_map(mod(ord-1, num_colors) + 1, :)); 
    hold on;
    if(settings.pH_status == 1) % pH decrescita lineare
        if(ord == 1)
            legend_str{ord} = sprintf('t = %d min, pH = %.2f', ord-1, pH(1));
        else
            legend_str{ord} = sprintf('t = %d min, pH = %.2f', ord-1, pH((ord-1)* Nstim));
        end
    elseif(settings.pH_status == 0) % pH fisso
        legend_str{ord} = sprintf('t = %d min, pH = %.2f', ord-1, pH(1));
    end
    legend(legend_str);
    grid on
    title(sprintf('APs shape'));
    xlim([0,settings.BCL+20]); % shift di 20 per visualizzazione ottimale
    xlabel('Time [ms]'); % il tempo qui è quello del PA, coincide con BCL (uguale per tutti)
    ylabel(sprintf('%s [%s]', names_sv(chosen), units_sv(chosen)));
    set(gca, 'FontSize',12,'LineWidth',1.5)

    % Capture the plot as a frame in the video
    frame = getframe(gcf);
    writeVideo(v, frame);

    % Reimposta i vettori per il prossimo ord
    V_div = [];
    time_div = [];

    % Aggiungi una pausa di 0.5 secondo tra ogni plot
    pause(0.5);
end

% Close the video writer
close(v);

clc;

%% con colorbar + grande

chosen = 25;
V_div = [];
time_div = [];
pH_segno = [];

BCL = settings.BCL / 1000; %[s]
Nstim = 60000 / settings.BCL; % [–]

% In queste variabili vengono salvati i potenziali di interesse e i relativi tempi
V_div_cell = cell(1, settings.TSim / 60000 +1);  % Inizializza il cell array con dimensioni corrispondenti al numero di iterazioni
time_div_cell = cell(1, settings.TSim / 60000 +1);  % Inizializza il cell array con dimensioni corrispondenti al numero di iterazioni

% Initialize the video writer
videoFileName = 'plot_video.mp4';
v = VideoWriter(videoFileName,  'MPEG-4');
v.FrameRate = 2; % Adjust frame rate as needed
v.Quality = 100;
open(v);

figure('Color', 'w', 'Position', get(0, 'Screensize')); % Full screen figure
num_colors = max(1, (settings.TSim / 60000 +1));
color_map_time = cool(num_colors);

% Assuming the pH values are provided in the 'pH' array
min_pH = min(pH);
max_pH = max(pH);
num_pH_colors = max(1, (settings.TSim / 60000 +1));
color_map_pH = flipud(cool(num_pH_colors));

for ord = 1:(settings.TSim / 60000 +1)
    if (ord == 1)
        for indice = 1:length(time)
            if (0<= time(indice) && time(indice) < BCL)
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    elseif (ord == 2)
        for indice = 1:length(time)
            if ((BCL < time(indice) && time(indice) <= BCL*Nstim))
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    elseif (ord > 2)
        for indice = 1:length(time)
            if ((ord-2)*BCL*Nstim < time(indice) && time(indice) <= (ord-1)*BCL*Nstim)
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    end

    % Aggiungi V_div al cell array
    V_div_cell{ord} = V_div;
    time_div_cell{ord} = time_div;

    % if(settings.pH_status == 1) % pH decrescita lineare
    %     if(ord == 1)
    %         pH_segno = [pH_segno, pH(1)];
    %     else
    %         pH_segno = [pH_segno, pH((ord-1)* Nstim)];
    %     end
    % elseif(settings.pH_status == 0) % pH fisso
    %         pH_segno = [pH_segno, pH(1)];
    % end

    % Map the current pH value to a color
    if(settings.pH_status == 1) % pH decrescita lineare
        current_pH = pH(min(ord, length(pH)));
    elseif(settings.pH_status == 0) % pH fisso
        current_pH = pH(1);
    end

    % Find the index for the current pH color
    pH_color_index = round((current_pH - min_pH) / (max_pH - min_pH) * (num_pH_colors - 1)) + 1;
    pH_color = color_map_pH(pH_color_index, :);
    

    % Plot una volta per ogni ord
    plot(1000*(time_div-time_div(1)), V_div, 'LineWidth', 4, 'Color', color_map_time(mod(ord-1, num_colors) + 1, :)); 
    hold on;

    grid on
    title(sprintf('APs shape'));
    xlim([0, settings.BCL+20]); % shift di 20 per visualizzazione ottimale
    xlabel('Time [ms]'); % il tempo qui è quello del PA, coincide con BCL (uguale per tutti)
    ylabel(sprintf('%s [%s]', names_sv(chosen), units_sv(chosen)));
    set(gca, 'FontSize', 30, 'LineWidth', 1.5)

    % Reimposta i vettori per il prossimo ord
    V_div = [];
    time_div = [];

    % Aggiungi una pausa di 0.5 secondo tra ogni plot
    pause(0.5);
end

% % Add the colorbar for time
% colormap(color_map_time);
% colorbar('Ticks', linspace(0, 1, num_colors), 'TickLabels', num2cell(0:num_colors-1));
% caxis([0, num_colors-1]);
% title(colorbar, 'Time step');
% set(gca, 'FontSize', 12, 'LineWidth', 1.5)


% Add the colorbar for pH
colormap(color_map_pH);
colorbar('Ticks', linspace(0, 1, num_pH_colors), 'TickLabels', (linspace(max_pH, min_pH, num_pH_colors)));
clim([6.2, 7.2]);
title(colorbar, 'pH');



% Close the video writer
close(v);

clc;

%%

figure('color','w')
plot(time_div_cell{1}/0.8,V_div_cell{1},'Color','k','LineWidth',4)
ylim([-100,50])
xlim([-0.2,+1.2])
set(gca,'xtick',[])


%% skip this section (AP alternativi)

chosen = 25;
V_div = [];
time_div = [];

BCL = settings.BCL / 1000; %[s]
Nstim = 60000 / settings.BCL; % [–]

% In queste variabili vengono salvati i potenziali di interesse e i relativi tempi
V_div_cell = cell(1, settings.TSim / 60000 +1);  % Inizializza il cell array con dimensioni corrispondenti al numero di iterazioni
time_div_cell = cell(1, settings.TSim / 60000 +1);  % Inizializza il cell array con dimensioni corrispondenti al numero di iterazioni


figure;

for ord = 1:(settings.TSim / 60000 +1)
    if (ord == 1)
        for indice = 1:length(time)

            if (0<= time(indice) && time(indice) < BCL)
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    elseif (ord == 2)
        for indice = 1:length(time)
            if ((BCL < time(indice) && time(indice) <= BCL*Nstim))
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    elseif (ord > 2)
        for indice = 1:length(time)
            if ((ord-2)*BCL*Nstim < time(indice) && time(indice) <= (ord-1)*BCL*Nstim)
                V_div = [V_div; StateControl(indice,25)];
                time_div = [time_div; time(indice)];
            end
        end
    end


    % Aggiungi V_div al cell array
    V_div_cell{ord} = V_div;
    time_div_cell{ord} = time_div;

    % Plot una volta per ogni ord
    plot(1000*(time_div-time_div(1)),V_div, 'LineWidth', 2.5, 'Color',  [0.2, 1, 0.4]); 
    hold on;
    
    grid on
    title(sprintf('SV%d: %s', chosen, names_sv(chosen)));
    xlim([0,settings.BCL+20]); % shift di 20 per visualizzazione ottimale
    xlabel('Time [ms]'); % il tempo qui è quello del PA, coincide con BCL (uguale per tutti)
    ylabel(sprintf('%s [%s]', names_sv(chosen), units_sv(chosen)));
    set(gca, 'FontSize',17,'LineWidth',2)

    % Reimposta i vettori per il prossimo ord
    V_div = [];
    time_div = [];

    % Aggiungi una pausa di 0.5 secondo tra ogni plot
    pause(0.5);
end

clc;

%% METRICs COMPUTATION

V = []; 

% metriche di interesse / utili
Vrest = [];
Vmax = [];
Vpp = [];
max_dV_dt = [];
APD90 = [];
APD50 = [];
t00 = [];

for col = 1:size(V_div_cell,2)
    
    % itero su V e ti, che cambiano a ogni iterazione
    V = V_div_cell{1,col};
    ti = (time_div_cell{1,col} - time_div_cell{1,col}(1,1))* 1000;

    % salvo a ogni iter le metriche in un vettore ad hoc
    
    % Vrest
    Vrest = [Vrest, V(1)]; % [mV]

    % Vmax
    Vmax = [Vmax, max(V)]; % [mV]

    % Vpp !!! usare Vmin
    Vpp = [Vpp, max(V) - V(1)]; % [mV]

    % max_dV_dt
    dV_dt = diff(V) ./ diff(ti); % [mV/ms] = [V/s] --> questo non lo salvo a ogni iter, non mi serve, salvo solo il max
    max_dV_dt = [max_dV_dt, max(abs(dV_dt))]; % [V/s]

    % ADP90
    V90 = max(V) - 0.9 * (max(V) - V(1)); % [mV] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    [~, ti_max] = max(V); % Calcolo dell'indice del massimo di V
    t0 = ti(dV_dt == max(abs(dV_dt(1:ti_max)))); % [ms] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    t90 = interp1(V(ti_max:end), ti(ti_max:end), V90, "nearest"); % [ms] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    APD90 = [APD90, t90 - t0]; % [ms]

    % ADP50
    V50 = max(V) - 0.5 * (max(V) - V(1)); % [mV] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    [~, ti_max] = max(V); % Calcolo dell'indice del massimo di V
    t0 = ti(dV_dt == max(abs(dV_dt(1:ti_max)))); % [ms] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    t50 = interp1(V(ti_max:end), ti(ti_max:end), V50, "nearest"); % [ms] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    APD50 = [APD50, t50 - t0]; % [ms]

    t00 = [t00, t0];

end 

% TRIANGULATION
TRI = APD90 - APD50;

% CI Coupling Interval
CI = settings.BCL; %[ms]

% DI Diastolic Interval
DI = CI-APD90; %[ms]

% Visualizza risultati
% Definisci il vettore tempo in minuti (da 0 a N)
tempo = (0:settings.TSim/60000);

% Crea una nuova figura
figure('Position',[100, 100, 800, 500],'Color','w');

% Plot di Vrest
subplot(2, 2, 1); % 2 righe, 2 colonne, primo subplot
plot(tempo, Vrest,'k', 'LineWidth', 4);
title('V_{rest}');
xlabel('Time [min]');
ylabel('V_{rest} [mV]');
grid on;
set(gca, 'FontSize',17,'LineWidth',2);

% Plot di Vpp
subplot(2, 2, 2); % 2 righe, 2 colonne, secondo subplot
plot(tempo, Vpp,'g', 'LineWidth', 4);
title('V_{pp}');
xlabel('Time [min]');
ylabel('V_{pp} [mV]');
grid on;
set(gca, 'FontSize',17,'LineWidth',2);

% Plot di max_dV_dt
subplot(2, 2, 3); % 2 righe, 2 colonne, terzo subplot
plot(tempo, max_dV_dt,'b', 'LineWidth', 4);
title('dV/dt_{max}');
xlabel('Time [min]');
ylabel('dV/dt_{max} [V/s]');
grid on;
set(gca, 'FontSize',17,'LineWidth',2);

% Plot di APD90
subplot(2, 2, 4); % 2 righe, 2 colonne, quarto subplot
plot(tempo, APD90,'r', 'LineWidth', 4);
title('APD90');
xlabel('Time [min]');
ylabel('APD90 [ms]');
grid on;
set(gca, 'FontSize',17,'LineWidth',2);

% Imposta il titolo generale per la figura
%sgtitle('Metrics Over Time');

%%
figure()
plot(DI, APD90, 'b', 'LineWidth', 2.5);
title('Restitution Curve');
xlabel('DI [ms]');
ylabel('APD_{90} [ms]');
grid on;
set(gca, 'FontSize',17,'LineWidth',2.5);

%% confronto indici

figure_size = [100, 100, 800, 900]; % [left, bottom, width, height]

% Create the figure with the specified size
figure('Position', figure_size);
subplot(2,1,1)
plot(time_div_cell{1}*1000, V_div_cell{1}, 'LineWidth', 2.5);
title('Steady-State AP Shape @ t = 0 min, pH = 7.2');
xlabel('time [ms]');
ylabel('V [mV]');
xlim([-100,800])
ylim([-80,50])
grid on;
set(gca, 'FontSize',17,'LineWidth',2);
subplot(2,1,2)
plot((time_div_cell{16}-time_div_cell{16}(1,1))*1000, V_div_cell{16}, 'LineWidth', 2.5);
title('Steady-State AP Shape @ t = 15 min, pH = 6.2');
xlabel('time [ms]');
ylabel('V [ms]');
xlim([-100,800])
ylim([-80,50])
grid on;
set(gca, 'FontSize',17,'LineWidth',2);
%%
% Create the figure with the specified size
figure();
plot(time_div_cell{1}*1000, V_div_cell{1}, 'LineWidth', 2.5, 'Color','k','LineStyle','--');
hold on
plot((time_div_cell{16}-time_div_cell{16}(1,1))*1000, V_div_cell{16}, 'LineWidth', 2.5,'Color','r');
xlabel('time [ms]');
ylabel('V [mV]');
xlim([-100,800])
ylim([-80,50])
grid on;
set(gca, 'FontSize',17,'LineWidth',2);
legend('t = 0 min, pH = 7.2', 't = 15 min, pH = 6.2')
%%

% Create the figure with the specified size
figure();
plot(time_div_cell{1}*1000, V_div_cell{1}, 'LineWidth', 5, 'Color','k','LineStyle','-');

ttttt = time_div_cell{1}*1000;
V1 = V_div_cell{1};


x1 = ttttt(1:length(ttttt)/2);
y1 = V1(1:length(ttttt)/2);

x2 = ttttt(length(ttttt)/2:end);
y2 = V1(length(ttttt)/2:end);

% Crea una nuova figura
figure;

% Plotta il primo segmento con la linea tratteggiata e il primo colore
plot(x1, y1, 'color', [0.2, 0.6, 0.8], 'LineWidth', 7); % Linea tratteggiata rossa
hold on;

% Plotta il secondo segmento con la linea tratteggiata e il secondo colore
plot(x2, y2, 'color', [0.8, 0.4, 0.2], 'LineWidth', 7); % Linea tratteggiata blu

xlabel('time [ms]');
ylabel('V [mV]');
xlim([-100,800])
ylim([-100,50])
set(gca, 'FontSize',17,'LineWidth',2);


%%

explawson = linspace(-75, -72, 16);

figure;
plot(tempo, Vrest, 'b','LineWidth', 2.5);
hold on
plot(tempo, explawson, 'ro', 'LineWidth', 2, 'MarkerFaceColor', 'r');
title('V_{rest}');
xlabel('Time [min]');
ylabel('V_{rest} [mV]');
grid on;
set(gca, 'FontSize',15,'LineWidth',2);
legend('Simulated with Maleckar','Experimental (adapted from Lawson, 2020)','best')

%%

figure('Color', 'w', 'Position', get(0, 'Screensize')); % Full screen figure
plot([0,100],[100,APD90(end)/APD90(1)*100],'r*-','LineWidth',5)
hold on
plot([0,100],[100,max_dV_dt(end)/max_dV_dt(1)*100],'b*-','LineWidth',5)
hold on
plot([0,100],[100,Vpp(end)/Vpp(1)*100],'g*-','LineWidth',5)
hold on
plot([0,100],[100,Vrest(end)/Vrest(1)*100],'k*-','LineWidth',5)
grid on
legend
title('Metrics Variation (%)');
xlim([-10,110])
xlabel('Acidosis [%]'); % il tempo qui è quello del PA, coincide con BCL (uguale per tutti)
legend('APD_{90}','dV/dt_{max}','V_{pp}','V_{rest}','location','southwest');
set(gca, 'FontSize', 30, 'LineWidth', 1.5)

%%

% Placeholder per BCL values
BCL_values = 700:200:1100; % BCL values from 500 ms to 1500 ms with 100 ms increments
APD90_values = [];
DI_values = [];

settings.TSim = 60000; % 1 min

settings.NumStim = floor(settings.TSim / settings.BCL);  % Number of action potentials to calculate
settings.storeLast = settings.NumStim; % I only want the very last one of each series of stimulation

settings.pH_status = 1; % pH decrescente
%settings.pH_status = 0; % pH fisso

settings.saving_status = 0; % salva ogni minuto
%settings.saving_status = 0; % salva tutto

% Loop through each BCL value
for BCL = BCL_values
    settings.BCL = BCL;
    [StateVars, Ti, Ns, pH] = maleckar_main(settings);

    % Supponiamo che 'Vm' sia il potenziale di membrana (da StateVars)
    Vm = StateVars(:, 25); % Indice 25 come esempio, potrebbe variare
    time = Ti;

    % Placeholder per APD90 e DI per questo BCL
    APD90 = [];
    DI = [];
    
    % Trova i picchi del potenziale d'azione per identificare i battiti
    [~, locs] = findpeaks(Vm, 'MinPeakHeight', -30); % supponiamo che il picco sia sopra -30 mV

    for i = 2:length(locs)
        % Estrai il potenziale d'azione del battito corrente
        AP_start = locs(i-1);
        AP_end = locs(i);
        V_AP = Vm(AP_start:AP_end);
        t_AP = time(AP_start:AP_end);

        % Calcola APD90
        APD90_value = t_AP(find(V_AP <= (V_AP(1) - 0.9 * (V_AP(1) - V_AP(end))), 1)) - t_AP(1);
        APD90 = [APD90, APD90_value];

        % Calcola DI
        if i > 2
            DI_value = time(locs(i)) - (time(locs(i-1)) - APD90(end-1));
            DI = [DI, DI_value];
        end
    end

    % Aggiungi i valori calcolati alla lista generale
    APD90_values = [APD90_values, APD90(2:end)];
    DI_values = [DI_values, DI];
end

% Plot della curva di restituzione
figure;
plot(DI_values*1000, APD90_values*1000, 'o-');
xlabel('Diastolic Interval (DI) [ms]');
ylabel('APD90 [ms]');
title('Restitution Curve');
grid on;
%%

% Placeholder per BCL values
BCL_values = 500:200:1100; % BCL values from 700 ms to 1100 ms with 200 ms increments
APD90_values = [];
DI_values = [];

settings.TSim = 900000; % 1 min

settings.pH_status = 1; % pH decrescente
settings.saving_status = 0; % salva ogni minuto

% Loop through each BCL value
for BCL = BCL_values
    settings.BCL = BCL;
    settings.NumStim = floor(settings.TSim / settings.BCL);  % Number of action potentials to calculate
    settings.storeLast = settings.NumStim; % I only want the very last one of each series of stimulation

    [StateVars, Ti, Ns, pH] = maleckar_main(settings);

    % Supponiamo che 'Vm' sia il potenziale di membrana (da StateVars)
    Vm = StateVars(:, 25); % Indice 25 come esempio, potrebbe variare
    time = Ti;

    % Placeholder per APD90 e DI per questo BCL
    APD90 = [];
    DI = [];

    % Trova i picchi del potenziale d'azione per identificare i battiti
    [~, locs] = findpeaks(Vm, 'MinPeakHeight', -30); % supponiamo che il picco sia sopra -30 mV

    for i = 2:length(locs)
        % Estrai il potenziale d'azione del battito corrente
        AP_start = locs(i-1);
        AP_end = locs(i);
        V_AP = Vm(AP_start:AP_end);
        t_AP = time(AP_start:AP_end);

        % Calcola APD90
        APD90_threshold = V_AP(1) - 0.9 * (V_AP(1) - min(V_AP));

    %     V0 = Vm(1);
    %     dV_dt = diff(V_AP) ./ diff(t_AP); % [mV/ms] = [V/s] --> questo non lo salvo a ogni iter, non mi serve, salvo solo il max 
    %     max_dV_dt =  max(abs(dV_dt)); % [V/s]
    %     V90 = max(V_AP) - 0.9 * (max(V_AP) - V0); % [mV] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    %     [~, ti_max] = max(V_AP); % Calcolo dell'indice del massimo di V
    %     t0 = ti(dV_dt == max(abs(dV_dt(1:ti_max)))); % [ms] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    % t90 = interp1(V(ti_max:end), ti(ti_max:end), V90, "nearest"); % [ms] --> questo non lo salvo a ogni iter, non mi serve, salvo solo ADP90
    % APD90 = [APD90, t90 - t0]; % [ms]


        APD90_index = find(V_AP <= APD90_threshold, 1);
        if ~isempty(APD90_index)
            APD90_value = t_AP(APD90_index) - t_AP(1);
            APD90 = [APD90, APD90_value];

        APP_start = locs(i);
        APP_end = locs(i);
        V_APP = Vm(APP_start:APP_end);
        t_APP = time(APP_start:APP_end);

            % Calcola DI
            if i > 1
               % DI_value = t_AP(1) - (t_AP(find(Vm(locs(i-1):AP_start) <= APD90_threshold, 1, 'last') + locs(i-1) - 1)); %%% ERRORE
               DI_value = t_APP(1) - t_AP(APD90_index) ; %%% ERRORE

                DI = [DI, DI_value];
            end
        end
    end

    % Aggiungi i valori calcolati alla lista generale
    APD90_values = [APD90_values, APD90];
    DI_values = [DI_values, DI];
end
%%
% Plot della curva di restituzione
figure;
plot(DI_values, APD90_values, 'o-');
xlabel('Diastolic Interval (DI) [ms]');
ylabel('APD90 [ms]');
title('Restitution Curve');
grid on;

%%

a=[275, 340, 360, 390, 460, 530, 550, 680, 730, 870,900,950,1000]
b=[110, 120, 130, 145, 150,170, 190, 200, 210, 220,225,227,227]

bs=[113, 135, 155, 170, 190,220, 250, 270, 280, 307,310,315,315]

figure('Color', 'w', 'Position', get(0, 'Screensize')); % Full screen figure
pp=interp1(a,b,200:10:1000,'linear');
pps=interp1(a,bs,200:10:1000,'linear');
plot(200:10:1000,pp,'o-','linewidth',5,'Color',  [1, 0.6, 0.6]);
hold on
plot(200:10:1000,pps,'o-','linewidth',5,'Color',[0.2, 1, 0.4]);
grid on
xlabel('Diastolic Interval (DI) [ms]');
ylabel('APD90 [ms]');
xlim([170,1030])
title('Restitution Curve');
legend('Acidotic','Healthy','location','northwest');
set(gca, 'FontSize', 30, 'LineWidth', 1.5)