%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                       %
%   In this M-file one can introduce new parameters into the model,     %
%   also diffrent initial conditions and times characterictic for       %
%   simulation may be chosen.                                           %
%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%############### Parametrization for system of ODE's

%############## Simulation time points ##################################


t0=0;                     % time when IFN is being introduced into the system
koniec=24;               % number of hours for simulation
tw=3600*koniec-1;         % length of IFN stimulation
te=3600*koniec;           % length of simulation from introducing IFN 


% Switching parameters
%Activation=0;
NTB=1;              % NTB=1 = no translation blocking 
                    % NTB=0 - translation blocking (CHX)
PB=0;               % PB=0 - phosphorylation blocking (SOCS1 or similar)
                    % PB=1 - no phosphorylation blocking
A_R=0;              % A_R=0 - no active repressor of IRF1 transcription
                    % A_R=1 - active repressor of IRF1 transcription
IA=0;               % IA=1 - inactive STAT1

block=1;             % if 0 then there is no saturation in phosphorylation; 1 otherwise

output=0;            % output=1 - concentrations, output=0 - number of molecules


% basic constants and rates
%---------------------------

Av=6.022e23;                  %Avogadro number

n_vol = 1.2e-12;              % nuclear volume HeLa cell - similar to Gorlich at al 2002
c_vol = n_vol/2;              % cytoplasm volume HeLa cell - similar to Gorlich at al 2002, corrected according to pictures

kv=c_vol/n_vol;

n_genes=4;                          % number of gene copies
v_max=n_genes*.16/Av/c_vol*1e6;     % [uM/s] - maximum transcription rate (per cytoplasmic volume)
v_transcription=.5e-5;
v_stat1_transcription=.6e-5;
v_stat2_transcription=0;%.1e-5;
v_lmp2_transcription = .2e-5;
v_tap1_transcription = .2e-5;
ktransl = 0.5;                      % [1/s] - translation rate for all transcripts

%ratio=0;          %percentage of STAT1 made inactive in the nucleus


% plot scaling
%---------------------------
if (koniec > 2)
    time_scale=3600;
else
    time_scale=60;
end

if (output == 0)
    cyt_scale=Av*c_vol/1e6;
    nuc_scale=Av*n_vol/1e6;
else
    cyt_scale=1;
    nuc_scale=1;
end

% time constants for inertial elements - in seconds

time_constant1=3600*.85;
time_constant2=3600*.85;
time_constant3=3600*.85;
time_constant4=3600*.85;
time_constant5=3600*.85;

time_constant_lmp=5*60;


constant1=1/time_constant1;
constant2=1/time_constant2;
constant3=1/time_constant3;
constant4=1/time_constant4;
constant5=1/time_constant5;
constant_lmp=1/time_constant_lmp;
% t4 = 3600*1.5;               % time constant
% t_2= t4*t4;
% t_3= t_2*t4;
% t_4= t_3*t4;

% Half-times of proteins and complexes [min]:
%--------------------------------------
               % transcripts
ts1t=.5*60;     % STAT1 mRNA
ts2t=.5*60;     % STAT2 mRNA
tirf1t=3*60;    % IRF1 mRNA
tlmp2t = 3*60;  % LMP2 mRNA
ttap1t = 2*60;  % TAP1 mRNA

               % proteins
ts1=15*60;     % STAT1
ts2=15*60;     % STAT2
ti1_in=20;   % inactive IRF1
ti1=20;        % active IRF1

ts1pc=20;      % STAT1p in cytoplasm (for dephosphorylation)
ts2pc=20;      % STAT2p in cytoplasm (for dephosphorylation)

               % complexes
ts1s1p_c=25;   % STAT1p|STAT1p complex in cytoplasm
ts1s1p_n=25;   % STAT1p|STAT1p complex in nucleus
ts1s2p_c=25;   % STAT1p|STAT2p complex in cytoplasm
ts1s2p_n=25;   % STAT1p|STAT2p complex in nucleus
tpiass1s1=3;   % PIAS|STAT1p|STAT1p
ts1i1 = 4;     % STAT1|IRF1 
ts1i1cyt = 30; % STAT1|IRF1 cytoplasmic


% degradation rates for single molecules [1/s]
%-----------------------------------------------------
                               % proteins
ks1deg = log(2)/ts1/60;        % STAT1 
ks1pdeg = ks1deg;           % STAT1p
ks2deg = log(2)/ts2/60;        % STAT2
ks2pdeg = ks1pdeg;          % STAT2p
ki1deg = log(2)/ti1/60;        % IRF1 active
ki1_indeg = log(2)/ti1_in/60;        % IRF1 active

                               % mRNA
kdegs1t=log(2)/ts1t/60;        % STAT1 mRNA 
kdegs2t=log(2)/ts2t/60;        % STAT2 mRNA 
kdegi1t=log(2)/tirf1t/60;      % IRF1t
kdeglmp2t = log(2)/tlmp2t/60;      % LMP2 mRNA
kdegtap1t = log(2)/ttap1t/60;      % TAP1 mRNA


% degradation rates for complexes (into substrates) [1/s]
%--------------------------------------------------------------------

% for dimers of phosphorylated STATs - it represents both dephosphorylation and degradation
kinvs1s1 = log(2)/ts1s1p_c/60;      % STAT1p|STAT1p
kinvs1s1n= log(2)/ts1s1p_n/60;      % (STAT1p|STAT1p)n
kinvs1s2 = log(2)/ts1s2p_c/60;      % STAT1p|STAT2p 
kinvs1s2n= log(2)/ts1s2p_n/60;      % (STAT1p|STAT2p)n

kinvpiass1s1 = log(2)/tpiass1s1/60; % PIAS|STAT1p|STAT1p 
kinvs1i1 = log(2)/ts1i1/60;         % STAT1|IRF1
kinvs1i1cyt = log(2)/ts1i1cyt/60;   % STAT1|IRF1 cytoplasmic

% degradation rates for complexes (non-reversible degradation) [1/s]
%--------------------------------------------------------------------

ks1s1pdeg = 0;0.05*kinvs1s1;                % (STAT1p|STAT1p) and (STAT1p|STAT1p)n
ks1s2pdeg = 0;0.1*kinvs1s2;                % (STAT1p|STAT2p) and (STAT1p|STAT2p)n
ks1i1deg = 0.01*kinvs1i1;                % STAT1|IRF1
ks1i1cytdeg = 0.01*kinvs1i1cyt;                % STAT1|IRF1 cytoplasmic
kdegpiass1s1 = 0;

% complexes creation rates [1/(uM*s)]
%---------------------------------------------
ks1s1=4;%log(2)/tc_s1s1/60;%.6;           % STAT1p|STAT1p 
ks1s2=10;%log(2)/tc_s1s2/60;%.9;          % STAT1p|STAT2p 
kpiass1s1 = .7;                           % PIAS|STAT1p|STAT1p
ks1i1 = .01;                                % STAT1|IRF1
ks1i1cyt = .001                             % STAT1|IRF1 cytoplasmic


%Activation of proteins
%--------------------------------------------------------------
kactivation = .7e-3;   % PIAS
tpias_inact = 20*60;
kpias_inact = log(2)/tpias_inact/60;

tacti1=1;
kacti1= log(2)/tacti1/60;
tinacti1=.25*60;
kinacti1=log(2)/tinacti1/60;


% Transcription thresholds [uM]
%---------------------------
kI1=6e-2;       % for IFN - induced IRF1 transcription
ks1=1e-3;



%Initial number of molecules
%corresponding to steady state in unstimulated cells
%-----------------------------------------------------

stat1_molecules = 1e5;                % number of STAT1 molecules in an unstimulated cell
stat2_molecules = stat1_molecules/1.5;    % number of STAT2 molecules in an unstimulated cell
pias_molecules = .5e5;                 % number of free inactive PIAS molecules in an unstimulated cell

lmpt2_molecules = 8;
tap1t_molecules = 8;

perc_s1_c = .9;                       % percentage of STAT1 molecules in cytoplasm - Assumption - 90%
perc_s2_c = .6;                      % percentage of STAT2 molecules in cytoplasm - Assumption - 99%


%Initial cellular concentrations - all [uM] (initial conditions are below them)
%corresponding to steady state in unstimulated cells
%-----------------------------------------------------

pias_ss = pias_molecules/n_vol/Av*1e6;
stat1c_ss = perc_s1_c*stat1_molecules/c_vol/Av*1e6;            % STAT1 in cytoplasm -  of molecules in cytoplasm
tot_stat1n_ss = (1-perc_s1_c)*stat1_molecules/n_vol/Av*1e6;    % total STAT1 in nucleus (including complexes)-

stat1n_ss = tot_stat1n_ss;                                     % free STAT1 in nucleus (separate variable
                                                               % since some complexes can be added in the future)

stat2c_ss = perc_s2_c*stat2_molecules/c_vol/Av*1e6;            % STAT2 in cytoplasm
stat2n_ss = (1-perc_s2_c)*stat2_molecules/n_vol/Av*1e6;        % STAT2 in nucleus

stat1t_ss = ks1deg/ktransl*(stat1c_ss  + stat1n_ss/kv);
                                   % STAT1 mRNA
                                   
stat2t_ss = ks2deg/ktransl*(stat2c_ss  + stat2n_ss/kv);

lmp2t_ss = lmpt2_molecules/c_vol/Av*1e6;
tap1t_ss = tap1t_molecules/c_vol/Av*1e6;
                                   
% Constitutive production rates [uM/s]
%------------------------------------
ks1tprod = kdegs1t*stat1t_ss;         % STAT1 mRNA
ks2tprod = kdegs2t*stat2t_ss;         % STAT2 mRNA
lmp2tprod = kdeglmp2t*lmp2t_ss;       % LMP2 mRNA
tap1tprod = kdegtap1t*tap1t_ss;       % TAP1 mRNA

% x protein is assumed to be produced only in stimulated cells
% PIAS is not degraded (actually assumed in the steady state concerning production and degradation),
% hence no constitutive production is needed

% Nuclear import/export rates [1/s]
%------------------------------------
transport=8;                   % [min]
es1=log(2)/transport/60;       % nuclear export of unphosphorylated STAT1 (Yamada)
kes1=0;
is1=stat1n_ss/stat1c_ss*(es1 + ks1deg/kv);
                               % nuclear import of unphosphorylated STAT1
                              
es2 = es1/1;                     % nuclear export of STAT2
is2 = (es2 + ks2deg/kv)*stat2n_ss/stat2c_ss;
                               % nuclear import of STAT2

transports1s1=.5;                % [min]
is1s1=log(2)/transports1s1/60;  % nuclear import of STAT1p|STAT1p
es1s1=0;               % nuclear export of STAT1p|STAT1p
transports1s2=transports1s1;    % it is done by the same importin alpha5
is1s2=log(2)/transports1s2/60;  % nuclear import of STAT1p|STAT2p
es1s2=0;               % nuclear export of STAT1p|STAT2p

transporti1=.5;
ii1 = log(2)/transporti1/60;    %nuclear import of active IRF1
ei1 = 1e-4*ii1;

transporti1_in=2;
ei1_in = log(2)/transporti1_in/60;    %nuclear export of inactive IRF1


%############### Initial conditions in [uM] ######################

y0=zeros(1,35);
y0(1)  = stat1c_ss;              % STAT1 unphosphorylated in cytoplasm
y0(2)  = stat2c_ss;              % STAT2 unphosphorylated in cytoplasm
y0(3)  = 0;                      % STAT1 phosphorylated in cytoplasm
y0(4)  = 0;                      % STAT2 phosphorylated in cytoplasm
y0(5)  = stat1n_ss;              % STAT1 unphosphorylated in nucleus 
y0(6)  = stat2n_ss;              % STAT2 unphosphorylated in nucleus 
y0(7) = 0;                      % (STAT1p|STAT1p) in cytoplasm
y0(8) = 0;                      % (STAT1p|STAT2p) in cytoplasm
y0(9) = 0;                      % (STAT1p|STAT1p) in nucleus
y0(10) = 0;                      % (STAT1p|STAT2p) in nucleus
y0(11) = stat1t_ss;              % STAT1 transcript     
y0(12) = 0;                      % active free PIAS
y0(13) = 0;                      % (PIAS|STAT1|STAT1)
y0(14) = pias_ss;                % inactive PIAS
y0(15)=0;                        % IRF1 mRNA
y0(16)=0;                        % IRF1 inactive protein
y0(17)=0;                        % IRF1 active in cytoplasm
y0(18)=0;                        % IRF1 active in nucleus
y0(19)=0;                        % IRF1 inactive in nucleus
y0(20)=0;
y0(21)=0;
y0(22)=0;
y0(23)=0;
y0(24)=stat2t_ss;
y0(25)=0;
y0(30)=lmp2t_ss;
y0(31)=tap1t_ss;
y0(32)=0;
y0(33)=0;
y0(34)=0;
y0(35)=0;


%Phosphorylation and dephosphorylation rates [1/s] 
%-------------------------------------------
% there is no dephosphorylation of STATs in nucleus due to assumption that there are no phosphorylated monomers there
% for a rate for dimers see the section with degradation of complexes into substrates
% phosphorylation, in turn, takes place only in cytoplasm
tphos1 = 32;           % phosphorylation of STAT1
tphos2 = 30;      % phosphorylation of STAT2

ks1phos  = log(2)/tphos1/60;           	  % phosphorylation rate for STAT1 - Yamada (?)
ks1dephc = log(2)/ts1pc/60;               % dephosphorylation of STAT1 in cytoplasm  - .5e-3 in Yamada (?)
ks2phos  = log(2)/tphos2/60; 	          % phosphorylation rate for STAT2
ks2dephc = log(2)/ts1pc/60;               % dephosphorylation of STAT2 in cytoplasm

ks1phos_sat = 1e4;                   %saturation of phosphorylation due to competition for receptors
ks2phos_sat = 1e4;

ks1phos=ks1phos*(1+block*ks1phos_sat*y0(1));
ks2phos=ks2phos*(1+block*ks2phos_sat*y0(2));