%comparison of simulation and data

%Experimental data notation
% time          - time points at which measures were taken
% stat1_n_data  - total STAT1 protein nuclear concentration
% stat1c_data   - total STAT1 protein cytoplasmic concentration
% stat1pc_data  - phosphorylated STAT1 in cytoplasm
% stat1pn_data  - phosphorylated STAT1 in nucleus
% time_hom      - time points at which EMSA measures of STAT1 homodimers were taken
% s1_s1_data    - Homodimer of phosphorylated STAT1
% irf1n_data    - Total IRF1 in nucleus
% irf1c_data    - Total IRF1 in cytoplasm

blot_data;  % reading the file with measurements
pcr_data;

% substracting blot background

%STAT1
stat1_n_data = stat1_n_data - min(stat1_n_data(length(stat1_n_data)), min(stat1_n_data));
stat1_n_data = stat1_n_data(1:length(stat1_n_data)-1);

stat1c_data = stat1c_data - min(stat1c_data(length(stat1c_data)), min(stat1c_data));
stat1c_data = stat1c_data(1:length(stat1c_data)-1);

stat1pc_data = stat1pc_data - min(stat1pc_data(length(stat1pc_data)), min(stat1pc_data));
stat1pc_data = stat1pc_data(1:length(stat1pc_data)-1);

stat1pn_data = stat1pn_data - min(stat1pn_data(length(stat1pn_data)), min(stat1pn_data));
stat1pn_data = stat1pn_data(1:length(stat1pn_data)-1);

% STAT2
stat2n_data = stat2n_data - min(stat2n_data(length(stat2n_data)), min(stat2n_data));
stat2n_data = stat2n_data(1:length(stat2n_data)-1);

stat2c_data = stat2c_data - min(stat2c_data(length(stat2c_data)), min(stat2c_data));
stat2c_data = stat2c_data(1:length(stat2c_data)-1);

stat2pc_data = stat2pc_data - min(stat2pc_data(length(stat2pc_data)), min(stat2pc_data));
stat2pc_data = stat2pc_data(1:length(stat2pc_data)-1);

stat2pn_data = stat2pn_data - min(stat2pn_data(length(stat2pn_data)), min(stat2pn_data));
stat2pn_data = stat2pn_data(1:length(stat2pn_data)-1);


s1_s1_data = s1_s1_data - min(s1_s1_data(length(s1_s1_data)), min(s1_s1_data));
s1_s1_data = s1_s1_data(1:length(s1_s1_data)-1);

s1_s2_data = s1_s2_data - min(s1_s2_data(length(s1_s2_data)), min(s1_s2_data));
s1_s2_data = s1_s2_data(1:length(s1_s2_data)-1);

irf1n_data = irf1n_data - irf1n_data(length(irf1n_data));
irf1n_data = irf1n_data(1:length(irf1n_data)-1);

irf1c_data = irf1c_data - irf1c_data(length(irf1c_data));
irf1c_data = irf1c_data(1:length(irf1c_data)-1);

isgf3_data = isgf3_data - isgf3_data(length(isgf3_data));
isgf3_data = isgf3_data(1:length(isgf3_data)-1);

time=time+5;
time(1)=0;
time_het=time_het+5;
time_het(1)=0;
stat2pc_time=stat2pc_time+5;
stat2pc_time(1)=0;

figure(10);
show_results(T,(Y(:,1)+Y(:,3)+2*Y(:,7)+Y(:,8))*cyt_scale,time,stat1c_data,koniec,3,2,1,'STAT1c');
show_results(T,(Y(:,3)+2*Y(:,7)+Y(:,8))*cyt_scale,time,stat1pc_data,koniec,3,2,2,'STAT1pc');
show_results(T,(Y(:,5)+2*Y(:,9)+Y(:,10)+2*Y(:,13)+Y(:,25))*nuc_scale,time,stat1_n_data,koniec,3,2,3,'STAT1n');
show_results(T,(2*Y(:,9)+Y(:,10)+2*Y(:,13))*nuc_scale,time,stat1pn_data,koniec,3,2,4,'STAT1pn');
show_results(T,Y(:,9)*nuc_scale,time_hom,s1_s1_data,koniec,3,2,5,'STAT1|STAT1');
if(koniec<=10)
	show_results(T,Y(:,10)*cyt_scale,time_het,s1_s2_data,koniec,3,2,6,'STAT1|STAT2');
end;
hold on;
if(koniec==10)
	plot(time_het/60, isgf3_data/trapz(time_het/60,isgf3_data), 'bx');
end
if(koniec==4)
   plot(time_het(1:7)/60, isgf3_data(1:7)/trapz(time_het(1:7)/60,isgf3_data(1:7)), 'bx');
end

figure(11);
show_results(T,(Y(:,2)+Y(:,4)+Y(:,8))*cyt_scale,time_het,stat2c_data,koniec,3,2,1,'STAT2c');
show_results(T,(Y(:,4)+Y(:,8))*cyt_scale,stat2pc_time,stat2pc_data,koniec,3,2,2,'STAT2pc');
show_results(T,(Y(:,6)+ Y(:,10))*nuc_scale,time_het,stat2n_data,koniec,3,2,3,'STAT2n');
show_results(T,(Y(:,10))*nuc_scale,time_het,stat2pn_data,koniec,3,2,4,'STAT2pn');
show_results(T,Y(:,10)*cyt_scale,time_het,s1_s2_data,koniec,3,2,6,'STAT1|STAT2');

figure(12);
subplot(3,3,1);
plot(T,Y(:,15)*cyt_scale);
grid on;
title('IRF1 mRNA (linear)');
show_results(T,Y(:,15)*cyt_scale,time_pcr,irf1t_data,koniec,3,3,2,'IRF1 mRNA');
subplot(3,3,4)
plot(T,Y(:,16)*cyt_scale);
grid on;
title('IRF1c inactive');
subplot(3,3,5)
plot(T,Y(:,17)*cyt_scale);
grid on;
title('IRF1c active');
subplot(3,3,3);
plot(T,Y(:,18)*nuc_scale);
grid on;
title('free IRF1n active')
subplot(3,3,6);
plot(T,Y(:,19)*nuc_scale);
grid on;
title('IRF1n inactive')

show_results(T,(Y(:,16)+Y(:,17))*cyt_scale,time,irf1c_data,koniec,3,3,7,'IRF1c total');
show_results(T,(Y(:,18)+Y(:,19)+Y(:,25))*nuc_scale,time,irf1n_data,koniec,3,3,8,'IRF1n total');

subplot(3,3,9)
plot(T, Y(:,25)*nuc_scale)
grid on;
title('STAT1|IRF1');

figure(13);
subplot(2,2,1);
plot(T,v_stat1_transcription*Y(:,23),'b-')
hold on;
plot([0 koniec], [v_max v_max],'r:')
grid on;
title('transcription rate');
legend('linear', 'max');

subplot(2,2,2);
plot(T,Y(:,11)*cyt_scale);
grid on;
title('STAT1 mRNA');

show_results(T,Y(:,11)*cyt_scale,time_pcr,stat1t_data,koniec,2,2,3,'STAT1 mRNA');

figure(14)
subplot(2,1,1);
plot(T, Y(:,30)*cyt_scale);
title('LMP2t');
show_results(T,Y(:,30)*cyt_scale,time_pcr,lmp2t_raw_data,koniec,2,1,2,'LMP2 mRNA');

% if (nuc_scale ~= 1)
%     figure(11);
%     plot(T, tot_S1_n+tot_S1_c);
%     title('Total STAT protein');
% end
% figure(12)
% plot(T, - ks1deg*Y(:,1) - ks1phos*Y(:,1) + ks1dephc*Y(:,3) - is1*Y(:,1) + es1*Y(:,5) + NTB*ktransl*Y(:,11) + 2*kinvs1s1*Y(:,7) + kinvs1s2*Y(:,8));
% hold on;
% plot(T, - ks1phos*Y(:,1) + ks1dephc*Y(:,3) - is1*Y(:,1) + es1*Y(:,5)  + 2*kinvs1s1*Y(:,7) + kinvs1s2*Y(:,8),'b--');
% title('Cytoplasmic, unphosphorylated STAT1 derivative')