close all

if (output == 0)
    label_y = 'Number of molecules';
else
    label_y = 'Molar concentration [uM]';
end

figure(1)
sim_scale=trapz(T, tot_S1_c);
data_scale1=trapz(time/60, stat1_c_g);
data_scale2=trapz(time/60, stat1c_data);

plot(T, tot_S1_c,'LineWidth',2);
ylabel('Number of molecules');
hold on;
plot(time/60, stat1_c_g/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time/60, stat1c_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT1 protein in cytoplasm (total)')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
sim_scale=trapz(T,tot_S1_pc);
data_scale1=trapz(time/60, stat1p_c_g);
data_scale2=trapz(time/60, stat1pc_data);

plot(T,tot_S1_pc,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time/60, stat1p_c_g/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time/60, stat1pc_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');
title('Phosphorylated STAT1 in cytoplasm')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3)
sim_scale=trapz(T,tot_S1_n);
data_scale1=trapz(time/60, stat1_n_g);
data_scale2=trapz(time/60, stat1_n_data);

plot(T,tot_S1_n,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time/60, stat1_n_g/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time/60, stat1_n_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT1 protein in nucleus (total)')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4)
sim_scale=trapz(T,tot_S1_pn);
data_scale1=trapz(time/60, stat1p_n_g);
data_scale2=trapz(time/60, stat1pn_data);

plot(T,tot_S1_pn,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time/60, stat1p_n_g/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time/60, stat1pn_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('Phosphorylated STAT1 in nucleus')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5)
sim_scale=trapz(T,Y(:,9)*nuc_scale);
data_scale2=trapz(time_hom/60, s1_s1_data);

plot(T,Y(:,9)*nuc_scale,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time_hom/60, s1_s1_data/data_scale2*sim_scale,'bo','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT1 homodimers in nucleus');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(6);
plot(T,Y(:,10)*nuc_scale);
title('STAT heterodimers in nucleus')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7)
% plot(T,(Y(:,2)+Y(:,4)+Y(:,8))*cyt_scale);
% title('STAT2 protein in cytoplasm (total)')

sim_scale=trapz(T, (Y(:,2)+Y(:,4)+Y(:,8))*cyt_scale);
data_scale2=trapz(time_het/60, stat2c_data);

plot(T, (Y(:,2)+Y(:,4)+Y(:,8))*cyt_scale,'LineWidth',2);
ylabel('Number of molecules');
hold on;
plot(time_het/60, stat2c_data/data_scale2*sim_scale,'bo','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT2 protein in cytoplasm (total)')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(8)
%plot(T, (Y(:,4)+Y(:,8))*cyt_scale);
%title('Phospho STAT2 in cytoplasm')
sim_scale=trapz(T, (Y(:,4)+Y(:,8))*cyt_scale);
data_scale2=trapz(stat2pc_time/60, stat2pc_data);

plot(T, (Y(:,4)+Y(:,8))*cyt_scale,'LineWidth',2);
ylabel('Number of molecules');
hold on;
plot(stat2pc_time/60, stat2pc_data/data_scale2*sim_scale,'bo','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT2 protein in cytoplasm (total)')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(9)
%plot(T, (Y(:,6)+ Y(:,10))*nuc_scale);
%title('STAT2 protein in nucleus (total)');
sim_scale=trapz(T, (Y(:,6)+Y(:,10))*nuc_scale);
data_scale2=trapz(time_het/60, stat2n_data);

plot(T, (Y(:,6)+Y(:,10))*nuc_scale,'LineWidth',2);
ylabel('Number of molecules');
hold on;
plot(time_het/60, stat2n_data/data_scale2*sim_scale,'bo','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT2 protein in nucleus (total)')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(10)
%plot(T,(Y(:,10))*nuc_scale);
%title('Phospho STAT2 in nucleus');
sim_scale=trapz(T, Y(:,10)*nuc_scale);
data_scale2=trapz(time_het/60, stat2pn_data);

plot(T, Y(:,10)*nuc_scale,'LineWidth',2);
ylabel('Number of molecules');
hold on;
plot(time_het/60, stat2pn_data/data_scale2*sim_scale,'bo','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('Phosporylated STAT2 in nucleus')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(11)
sim_scale=trapz(T,Y(:,15)*cyt_scale);
data_scale1=trapz(time_pcr/60, irf1t_raw_data);
data_scale2=trapz(time_pcr/60, irf1t_data);

plot(T,Y(:,15)*cyt_scale,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time_pcr/60, irf1t_raw_data/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time_pcr/60, irf1t_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('IRF1 mRNA');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(12)
sim_scale=trapz(T,(Y(:,16)+Y(:,17))*cyt_scale);
data_scale1=trapz(time/60, irf1_cyt_g);
data_scale2=trapz(time/60, irf1c_data);

plot(T,(Y(:,16)+Y(:,17))*cyt_scale,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time/60, irf1_cyt_g/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time/60, irf1c_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');

ylabel('Normalized results');
title('IRF1 protein in cytoplasm');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(13)
sim_scale=trapz(T,(Y(:,18)+Y(:,19)+Y(:,25))*nuc_scale);
data_scale1=trapz(time/60, irf1_nuc_g);
data_scale2=trapz(time/60, irf1n_data);

plot(T,(Y(:,18)+Y(:,19)+Y(:,25))*nuc_scale,'LineWidth',2)
ylabel(label_y);
hold on;
%plot(time/60, irf1_nuc_g/data_scale1*sim_scale,'bo','LineWidth',2);
plot(time/60, irf1n_data/data_scale2*sim_scale,'bo','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('IRF1 protein in nucleus');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(14)
sim_scale=trapz(T,Y(:,11)*cyt_scale);
data_scale1=trapz(time_pcr/60, stat1t_raw_data);
data_scale2=trapz(time_pcr/60, stat1t_data);

plot(T,Y(:,11)*cyt_scale,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time_pcr/60, stat1t_raw_data/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time_pcr/60, stat1t_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('STAT1 mRNA');
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(15)
sim_scale=trapz(T,Y(:,30)*cyt_scale);
data_scale1=trapz(time_pcr/60, lmp2t_raw_data);
data_scale2=trapz(time_pcr/60, lmp2t_data);

plot(T,Y(:,30)*cyt_scale,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time_pcr/60, lmp2t_raw_data/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time_pcr/60, lmp2t_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');

title('LMP2 mRNA');
grid on;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(16)
sim_scale=trapz(T,Y(:,31)*cyt_scale);
data_scale1=trapz(time_pcr/60, tap1t_raw_data);
data_scale2=trapz(time_pcr/60, tap1t_data);

plot(T,Y(:,31)*cyt_scale,'LineWidth',2);
ylabel(label_y);
hold on;
plot(time_pcr/60, tap1t_raw_data/data_scale1*sim_scale,'bo','LineWidth',2);
%plot(time_pcr/60, tap1t_data/data_scale2*sim_scale,'kx','LineWidth',2);

ax1=gca;
y1_limit=get(ax1,'YLim');
y1_limit=[0,y1_limit(2)];
set(ax1,'YLim',y1_limit,'XLim',[0,horizon],'YTick',[0,.25*y1_limit(2),.5*y1_limit(2),.75*y1_limit(2),y1_limit(2)]);
y2_limit=y1_limit/sim_scale;

axes('XAxisLocation','bottom','YAxisLocation','right',...
    'XLim',[0,horizon],'YLim',y2_limit,...
    'YTick',[0,.25*y2_limit(2),.5*y2_limit(2),.75*y2_limit(2),y2_limit(2)],...
    'Color','none');
ylabel('Normalized results');
title('TAP1 mRNA');
grid on;