clear all
close all
clc

%% variabili globali
global V1 V2 V3 t_dop y i ke3 k31 alpha beta gamma
global S Wgc Wgs Wnc Wns Ke STN_ON T_ON 
global c3_delay tfreq Tapping_f Dop_tonic 


%% inizializzazione sinapsi
Nc = 4;
load W_tot_new_W0e5_D1e0
    Wgc = squeeze(Wgc_epocs(:,:,100));
    Wgs = squeeze(Wgs_epocs(:,:,100));
    Wnc = squeeze(Wnc_epocs(:,:,100));
    Wns = squeeze(Wns_epocs(:,:,100));
Ke = 7;


%% carico i dati e i parametri del paziente
nome_paziente = input('Name of the patient data file (between apices)? ');
stringa = ['load ' nome_paziente];
nome_parametri = input('Name of the parameter file (between apices)? ');
stringa1 = ['load ' nome_parametri];
eval(stringa)
eval(stringa1)

%% calcolo la cinetica della levodopa
V1 = 12;
V2 = 32;
V3 = 2;
k31 = 0.02;
%ke3 = 0.03;

i = [zeros(Delay_levodopa,1);  3.33*(ones(300,1)); zeros(2200 - Delay_levodopa,1)]';

ke1 = ketot - k31;


dt = 0.1;
t1 = [0:dt:250];
L = length(t1);
c1 = zeros(L,1);   % plasma+periferico
c2 = zeros(L,1);   
c3 = zeros(L,1);
%c1(1) = D/V1;


for j = 1: L-1,
    dc1 = -k21/V1*c1(j)+k12/V1*c2(j)-ke1/V1*c1(j)-k31/V1*c1(j)+i(j)/V1;
    dc2= k21/V2*c1(j)-k12/V2*c2(j);
    c1(j+1) = c1(j) +dt*dc1;
    c2(j+1) = c2(j) + dc2*dt;
end

index = [1 151 301 451 601 751 901 1201 1501 1801];  
kk = length(y);
index = index(1:kk);  % inserito perché in qualche caso y è meno lungo

figure
width = 1.5;
font = 18;
plot(t1(index),c1(index),'bo-',t_dop,y,'r*--','linewidth', width)
xlabel('time (min)','fontsize',font)
ylabel('Levodopa concentration (\mu g/mL)','fontsize',font) 
set(gca,'fontsize',font)

%% parametri dei gangli della base
alpha = 0.75;  %(0.2*(Ugo_trigger-0.8)+0.5)/(0.7*(Ugo_trigger-0.8));

%guadagno da DA a No-Go (inibizione)
beta = -1;

%guadagno da DA a interneurone colinergico (inibizione)
gamma = -0.5;  

Ns = 4;
S = zeros(Ns,1);
S(1) = 1;
   

STN_ON = 1;
T_ON = 1;

dt = 0.1;
tau=30;
if ~exist('Nprove','var')
    Nprove=Np;
end

for k1 = 1: Nprove
        
Delay = p1_totale(k1);
ke3 = p2_totale(k1);
c3 = zeros(L,1);   

for j = 1: L-1,
    dc3= k31/V3*c1(j) -ke3/V3*c3(j);
    c3(j+1) = c3(j) + dc3*dt;
end

Delay_indice = round(Delay/dt);

% variabile c3_ritardata
c3_delay = [zeros(1,Delay_indice) c3(1:L-Delay_indice)']';
%il passo di integrazione per il modello della levodopo era 0.1. Per cui
% ricampiono a un minuto;
c3_delay = c3_delay(1:10:length(t1));

Dop_max = p3_totale(k1);
Dop_50 = p4_totale(k1);
N = p5_totale(k1);
Dop_tonic = p6_totale(k1);

DA_tot = [];
ft_tot = [];
for iiii=1:1:length(tfreq)
    index = tfreq(iiii)+1;
    Dop_ex = c3_delay(index);
    % assegno dopamina tot
    DA = Dop_tonic+(Dop_max*Dop_ex^N)/(Dop_50^N+Dop_ex^N);
    [Uc,C,Ugo,Go,IGo_DA_Ach,Unogo,NoGo,INoGo_DA_Ach,Ugpe,Gpe,Ugpi,Gpi,Ut,T,Ustn,STN,E,tt,k_tap_vett,Uchi,ChI,ft] = BG_model_function_tapping_mauro_3(S,Wgc,Wgs,Wnc,Wns,Ke,STN_ON,T_ON,DA);
    
    DA_tot = [DA_tot DA];
    ft_tot = [ft_tot ft];
    
end

%% faccio le figure
width = 1.5;
font = 18;



figure
plot(tfreq,ft_tot.*60,'bo-',tfreq,Tapping_f,'r*--','linewidth', width)
xlabel('time (min)','fontsize',font)
ylabel('tapping frequency (taps/min)','fontsize',font) 
set(gca,'fontsize',font)

    end