%         full model with IRF1|STAT1 complexes and LMP2 transcription,
    %         without time delay for STAT mRNA production (about 30 min)
    %         so data for STAT1 mRNA and partly STAT1 protein should appear
    %         a bit before actual data
    %         CBP proteins and IRF1|CBP complexes not modelled explicitely
    %         but included in the inertial element - up to 4th order
    %
    %
    %         y(1)   STAT1 unphosphorylated in cytoplasm
    %         y(2)   STAT2 unphosphorylated in cytoplasm
    %         y(3)   STAT1 phosphorylated in cytoplasm
    %         y(4)   STAT2 phosphorylated in cytoplasm
    %         y(5)   STAT1 unphosphorylated in nucleus
    %         y(6)   STAT2 unphosphorylated in nucleus
    %         y(7) (STAT1p|STAT1p) in cytoplasm
    %         y(8) (STAT1p|STAT2p) in cytoplasm
    %         y(9) (STAT1p|STAT1p) in nucleus
    %         y(10) (STAT1p|STAT2p) in nucleus
    %         y(11) STAT1 mRNA (linear)
    %         y(12) active, free PIAS
    %         y(13) PIAS complex with STAT1|STAT1
    %         y(14) inactive PIAS
    %         y(15) IRF1 mRNA
    %         y(16) IRF1 inactive in cytoplasm
    %         y(17) IRF1 active in cytoplasm
    %         y(18) IRF1 active in nucleus
    %         y(19) IRF1 inactive in nucleus
    %         y(20) inertial element (1st)
    %         y(21) inertial element (2nd)
    %         y(22) inertial element (3rd)
    %         y(23) inertial element (4th)
    %         y(24) STAT2 mRNA
    %         y(30) LMP2 mRNA
    %         y(31) TAP1 mRNA
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 function dy=model_pias(t,y,TR)

 parameters_pias; %parametrization included in M-file 'parameters.m'
 

%######################################################################### 

% TR=0 IFN off, TR=1 IFN on
dy=zeros(31,1);

% equation 1 - STAT1 unphosphorylated in cytoplasm
dy(1) = - ks1deg*y(1) - TR*ks1phos*y(1)/(1+block*ks1phos_sat*y(1)) + ks1dephc*y(3) - is1*y(1) + es1*y(5) + NTB*ktransl*y(11) + 2*kinvs1s1*y(7) + kinvs1s2*y(8);


% equation 2 - STAT2 unphosphorylated in cytoplasm
dy(2) = NTB*ktransl*y(24) - ks2deg*y(2) - TR*ks2phos*y(2)/(1+block*ks2phos_sat*y(2)) + ks2dephc*y(4) - is2*y(2) + es2*y(6) + kinvs1s2*y(8) ;


% equation 3 - STAT1 phosphorylated in cytoplasm
dy(3) = TR*ks1phos*y(1)/(1+block*ks1phos_sat*y(1)) - ks1pdeg*y(3) - ks1dephc*y(3) - 2*ks1s1*y(3)*y(3) - ks1s2*y(3)*y(4);


% equation 4 - STAT2 phosphorylated in cytoplasm
dy(4) = TR*ks2phos*y(2)/(1+block*ks2phos_sat*y(2)) - ks2pdeg*y(4) - ks2dephc*y(4) - ks1s2*y(3)*y(4);  



% equation 5 - STAT1 unphosphorylated in nucleus
dy(5) = is1*kv*y(1) - es1*kv*y(5) - ks1deg*y(5)  +  kinvs1s2n*y(10)  + 2*kinvs1s1n*y(9)+ 2*kinvpiass1s1*y(13) - ks1i1*y(18)*y(5) + kinvs1i1*y(25);


% equation 6 - STAT2 unphosphorylated in nucleus
dy(6) = is2*kv*y(2) - es2*kv*y(6) - ks2deg*y(6) + kinvs1s2n*y(10);


% equation 7 - STAT1p|STAT1p complex in cytoplasm; 
dy(7) = ks1s1*y(3)*y(3) - (ks1s1pdeg)*y(7) - is1s1*y(7) - kinvs1s1*y(7);



% equation 8 - STAT1p|STAT2p complex in cytoplasm; 
dy(8) = ks1s2*y(3)*y(4) - (ks1s2pdeg + is1s2)*y(8) - kinvs1s2*y(8);


% equation 9 - STAT1p|STAT1p complex in nucleus; 
dy(9) = - (kinvs1s1+ks1s1pdeg)*y(9) + is1s1*kv*y(7) - kpiass1s1*y(9)*y(12);


% equation 10 - STAT1p|STAT2p complex in nucleus; 
dy(10) = - (kinvs1s2+ks1s2pdeg)*y(10) + is1s2*kv*y(8);


% equation 11 - STAT1 transcript     
dy(11) = NTB*ks1tprod + v_stat1_transcription*y(23) - kdegs1t*y(11);


% equation 12 - free active PIAS
dy(12) = kactivation*y(10)*y(14) - kpiass1s1*y(9)*y(12) + kinvpiass1s1*y(13);


% equation 13 - PIAS|STAT1|STAT1
dy(13) = kpiass1s1*y(9)*y(12) - kinvpiass1s1*y(13) - kdegpiass1s1*y(13);


% equation 14 - inactive PIAS; 
dy(14) = -kactivation*y(10)*y(14);

% equation 15 - IRF1 transcript with transription rate linearly dependent on TF concentration
dy(15) = v_transcription*y(9) - kdegi1t*y(15);

% equation 16 - IRF1 inactive in cytoplasm
%dy(16) = NTB*ktransl*y(15) - ki1_indeg*y(16) - kacti1*y(16) + kinacti1*y(17);
dy(16) = - ki1_indeg*y(16) + ei1_in*y(19);

% equation 17 - IRF1 active in cytoplasm
%dy(17) = kacti1*y(16) - kinacti1*y(17) - ki1deg*y(17) - ii1*y(17) + ei1*y(18);
dy(17) = NTB*ktransl*y(15) - ki1deg*y(17) - ii1*y(17) + ei1*y(18);

% equation 18 - IRF1 active in nucleus
%dy(18) = kv*ii1*y(17) - kv*ei1*y(18) - ki1deg*y(18);
dy(18) = kv*ii1*y(17) - kv*ei1*y(18) - ki1deg*y(18) - kinacti1*y(18)  - ks1i1*y(18)*y(5) + kinvs1i1*y(25);

% equation 19 - IRF1 inactive in nucleus
dy(19) = kinacti1*y(18) - ki1_indeg*y(19) - kv*ei1_in*y(19);

% equation 20 - inertial element (first - with active IRF1n as an input)
dy(20)= - constant1*y(20) + constant1*y(18);

% equation 21 - inertial element (2nd)
dy(21)= - constant2*y(21) + constant2*y(20);

% equation 22 - inertial element (3rd)
dy(22)= - constant3*y(22) + constant3*y(21);

% equation 23 - inertial element (4th - last one - it is an input for STAT1 transcription)
dy(23)= - constant4*y(23) + constant4*y(22);

% equation 24 - STAT2 mRNA
dy(24)= NTB*ks2tprod + NTB*v_stat2_transcription*y(22) - kdegs2t*y(24);

% equation 25 - STAT1|IRF1 complexes in nucleus
dy(25) = ks1i1*y(18)*y(5) - (ks1i1deg + kinvs1i1)*y(25);

% equation 26 - inertial element for LMP2 transctription (first - with STAT1|IRF1 as an input)
dy(26)= - constant_lmp*y(26) + constant_lmp*y(25);

% equation 27 - inertial element (2nd)
dy(27)= - constant_lmp*y(27) + constant_lmp*y(26);

% equation 28 - inertial element (3rd)
dy(28)= - constant_lmp*y(28) + constant_lmp*y(27);

% equation 29 - inertial element (4th - last one - it is an input for LMP2 transcription)
dy(29)= - constant_lmp*y(29) + constant_lmp*y(28);

%equation 30 - LMP2 transcript
dy(30) = NTB*lmp2tprod + NTB*v_lmp2_transcription*y(25) - kdeglmp2t*y(30);

%equation 31 - TAP1 transcript
%inertial elements neglected - time constant quite low anyway,, as
%indicated by lmp2 results
dy(31) = NTB*tap1tprod + NTB*v_tap1_transcription*y(29) - kdegtap1t*y(31);