%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   11/29/04                                                    %
%   Code for solving system of ODE's given by function          %
%   model using 4th order matlab solver.                        %
%   Simulation  takes four steps.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;

parameters_pias;            % parametrization of ODE's and simulation times
tspan=[0,3600*koniec];
Activation=1;

% Step_one - before detachment of enhancewosome from STAT1 promoter site

[T,Y]=ode23tb(@model_pias,tspan,y0,[],Activation); 

%Step two - not much of STAT1 production
%----------------------------------
% 
% Activation = 1;
% tspan2=[3600*koniec1,3600*koniec];               
% x=size(Y);
% Y01=Y(x(1),:);         %transfer of initial condition from step one to step two
% [T1,Y1]=ode23tb(@model_1,tspan2,Y01,[],Activation); 
% 
% T=[T;T1];
% Y=[Y;Y1];


T=T/time_scale;

tot_S1_pn = (2*Y(:,9)+Y(:,10)+2*Y(:,13))*nuc_scale;
tot_S1_pc = (Y(:,3)+2*Y(:,7)+Y(:,8))*cyt_scale;
tot_S1_n = (Y(:,5)+2*Y(:,9)+Y(:,10)+2*Y(:,13)+Y(:,25))*nuc_scale;
tot_S1_c = (Y(:,1)+Y(:,3)+2*Y(:,7)+Y(:,8))*cyt_scale;

figure(1);
set(gcf,'Color',[1,1,1])

subplot(3,4,1);
plot(T,Y(:,1)*cyt_scale);
grid on;
title('STAT1');

subplot(3,4,2);
plot(T,Y(:,2)*cyt_scale);
grid on;
title('STAT2');

subplot(3,4,3);
plot(T,Y(:,3)*cyt_scale);
grid on;
title('STAT1p');


subplot(3,4,4);
plot(T,Y(:,4)*cyt_scale);
grid on;
title('STAT2p');


subplot(3,4,5);
plot(T,Y(:,7)*cyt_scale);
grid on;
title('1p|1p');

subplot(3,4,6);
plot(T,Y(:,8)*cyt_scale);
grid on;
title('1p|2p');

subplot(3,4,7);
plot(T,Y(:,9)*nuc_scale);
grid on;
title('(1p|1p)n');

subplot(3,4,8);
plot(T,Y(:,10)*nuc_scale);
grid on;
title('(1p|2p)n');

subplot(3,4,9);
plot(T,tot_S1_n);
grid on;
title('tot S1n');

subplot(3,4,10);
plot(T,tot_S1_c);
grid on;
title('tot S1c');

subplot(3,4,11);
plot(T,tot_S1_pn);
grid on;
title('tot S1pn');

subplot(3,4,12);
plot(T,tot_S1_pc);
grid on;
title('tot S1pc');

figure(3);

plot(T,v_transcription*Y(:,9),'g--');
hold on;
plot(T,v_max*Y(:,9)./(kI1 + Y(:,9)),'b');
hold on;
plot([0 koniec], [v_max v_max],'r:')
grid on;
title('transcription rate');
legend('linear', 'M-M', 'max');

compare;
% figure(6)
% 
% subplot(2,1,1);
% plot(T,(Y(:,12)+Y(:,13))*nuc_scale);
% grid on;
% title('PIAS active');
% 
% subplot(2,1,2);
% plot(T,Y(:,14)*nuc_scale);
% grid on;
% title('PIAS inactive');