% ODEs_t7vEB1 - ODEs for implementation of t7v1 from Endy1997(Biotech&Bioeng)
%
% Programmers:
%   Elsa Birch - ewbirch@stanford.edu
%   Nicholas Ruggero - nruggero@stanford.edu
%
% Documentation:
%   References to Endy1997B&B, Endy1997Thesis
%
% Dependancies:
%   Called by: setup_t7vEB1.m
%   Calls: 
%
% Chemical kinetic (ordinary differential equation model) of bacteriophage
% T7 replication in host E. coli
%
% Notes:
%   - Assumes temperature is 30C
%   - Time unit used is 1 second
%   - T7 genetic elements named with decimal converted to 'o' in variable
%   names, eg 1.1A becomes 1o1A


function [varargout] = T7_37Growth_EB_VirusODEs(t,C,varargin)
% t
% Standard varargin formatting setup for use with vFBA functions
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% varargin is expected as either non-existant or a cell array. In the case when the first
% cell of varargin is string 'identifiers' varargout is cell array of rxn and metabolite
% IDs corresponding to function outputs. Any other elements of varargin may be used by the
% function, and note that in the case of an identifiers request the cells of varargin are
% shifted after first cell is removed, so that with or without an identifier request the
% varargin format passed to the rest of the function is the same.


if (~isempty(varargin) && ischar(varargin{1})) && strcmpi( varargin{1}, 'identifiers' )
    % case when the first vararg in is identifier request flag
    produceIdentifiersFlag = true;
    varargin = varargin(2 : end);
else
    produceIdentifiersFlag = false; 
end % first cell of vararg in is an identifier request

% model is the first varargin, limits is the second
model = varargin{ 1 };
viralLimits = varargin{ 2 };

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if ~produceIdentifiersFlag
    
    
    % Prepare for rate calculations, book keeping only
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % This portion of the code only contains variable renamings or
    % initializations and no calculations from references
    
    % Divide up concentration arrays into mRNA conecentrations and gene product
    % (gp) concentrations so that they are more easily passed and referenced by index
    CmRNA = C((model.n_gene + 4 + 1):(2*model.n_gene + 4)); % Concentration of mRNA
    Cgp = C(1:model.n_gene); % Concentration of gene product (gp)
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ Prepare for rate calculations, book keeping only ]
    
    
    
    % Genome translocation - calculation of X
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % X(t) = length of genome inserted in host cell at time t
    X = translocation(t,model); % nt
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ genome translocation - calculation of X ]
    
    
    
    % Rate calculations for mRNA
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % Prepares all needed variables and evaluates ODEs for production of mRNA
    % species
    [dmRNAdt,~,mrnaProd] = dmRNAdtCalculator(t,X,Cgp,CmRNA,model,viralLimits);
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ Rate calculations for mRNA  ]
    
    % Rate limitation for phage assembly
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    
    % Rate Calculation for Progeny Phage
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    [cRateLimitingDNApackage, rateLimitPhageAssembly ] = determineT7assemblyRateLimit(C, model);
    
    % Calculates rate of phage formation as described
    % in Endy1997(Biotech&Bioeng) using Equation 10.
    %
    % dPhage
    % ------ = ( N_progeney ) * ( k_packagingT7DNA ) * P_R
    %   dt
    %
    % Here N_progeney is efficieny of DNA or Procapsids (P_R)
    %                 resulting in functional progeney, from ENdy1997B&B
    
    phageProd = model.Nprogeney * model.k_packagingT7DNA * cRateLimitingDNApackage;
    
    % if a component other than PC or DNA is limiting (phage stoich)
    % this assumes that the rate of phage assembly does not depend on
    % included protiens unless there are too few stoichiometrically
    % avalible in which case phage assembly rate is limited strictly to the
    % maximum possible with avalible components
    if phageProd > rateLimitPhageAssembly
        phageProd = rateLimitPhageAssembly;
    end
    
    % FBA Limits if virus draws more than is feasible
    if viralLimits{1}(3*model.n_gene+7) && (phageProd > viralLimits{2}(3*model.n_gene+7))
        phageProd = viralLimits{2}(3*model.n_gene+7);
    end % limit viral flux and over limit
    
    phageCons = 0;
    dPhagedt = phageProd - phageCons;
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    
    
    % Rate Calculations for Gene Products
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % Prepares all needed variables and evaluates ODEs for production, decay,
    % and use of gene products (gp)
    [dGPdt,gpProd] = dGPdtCalculator(Cgp,CmRNA,C,X,model,viralLimits, phageProd );
    % Also outputs, CrateLimitingDNAPackage, the concentration of the rate
    % limiting component of the T7 virus progeney (i.e. the component with the
    % lowest concentration of T7DNA and PC), which is used in later setps
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  [ Rate Calculations for Gene Products ]
    
    
    
    % Rate Calculations for Phage DNA
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % first calculate number of possible replisomes and replication forks
    % [You 2002]
    
    gp5 = Cgp(model.genomeT7.gp5.n_gene); % gp5 is the DNA polymerase
    gp4hex = Cgp(model.genomeT7.gp4a.n_gene) / 6; % gp4 functions as a hexamer helicase and primase
    replisomes = min( [gp4hex,gp5] ); % /cell
    DNA = C(model.n_gene+1);
    replicationForks = 2 * (DNA + 1);
    rep = min([replisomes,replicationForks]);
    
    
    % Binding of gp1 by gp 3.5 
    gp3o5boundT7RNAP = bindingEquilib( [Cgp(model.genomeT7.gp1.n_gene),Cgp(model.genomeT7.gp3o5.n_gene)],...
                        model.K_eqBindGp3o5, model.vol_Ecoli ,t);
    gp1_3o5 = gp3o5boundT7RNAP;
    gp2o5 = Cgp(model.genomeT7.gp2o5.n_gene);
    
    dNTPnonzero = ((C(model.n_gene + 4) / 2) > 0);
    thresholdYouGn = 5e-7 *6.022e23 * model.vol_Ecoli;
    replicateDNA = dNTPnonzero && (gp2o5 > 1) && (gp1_3o5 > thresholdYouGn); % also, is this in M? move constant
    
    % Note, that in Endy & You calcs are in bp/sec not nt/sec, however for
    % interaction with FBA this must be in terms of nt/sec (hence factor of 2)
    bpOfDNTPs = C(model.n_gene + 4) / 2;
    
    % Calculate the concentration rate of change of T7 DNA using Equation 7 from Endy1997(Biotech&Bioeng)
    %
    %  d(DNA)        [bp of dNTP] * k_repT7RNA * [rep]
    % ------- = -----------------------------------------------  - ( k_packagingT7DNA ) * P_R
    %   dt       ([bp of dNTP] + K_mDNAelongat) * (L_T7genome)
    %

    if replicateDNA
        % dna production will be rescaled so that resource/rxn is on same
        % order as other viral rxns
        dnaProd =  (bpOfDNTPs * model.k_replicationT7DNA * rep /...
            ( bpOfDNTPs  + model.K_mDNAelongation ) / model.lengthT7) * 40; 
    else
        dnaProd = 0; 
    end
    
    % FBA Limits if virus draws more than is feasible
    if viralLimits{1}(3*model.n_gene+5) && (dnaProd > viralLimits{2}(3*model.n_gene+5))
        dnaProd = viralLimits{2}(3*model.n_gene+5);
    end % limit viral flux and over limit
    
%     dnaCons = model.k_packagingT7DNA * CrateLimitingDNApackage;
    dnaCons = phageProd / model.Nprogeney;

    dDNAdt = dnaProd/40 - dnaCons;

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % - [Rate Calculations for Phage DNA]
    
    
% %     % Host genome degradation
% %     % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% %     
% %     % Calculates degradation of host DNA, as described but not explicitly stated
% %     % in [Endy1997(Biotech&Bioeng)]. Constant rate of production/degradation
% %     % k_degradationEcoliDNA.
% %     %
% %     % NOTE: Note the production time limits in model structure (only occurs
% %     % between 7.5 min and 15 min after transfection begins).
% %     %
% %     %                    k_degradationEcoliDNA        if 15min > t > 7.5min
% %     % dNTPproduction = {
% %     %                    0                            otherwise
% %     %
% %     hostGenomeDegrade = model.k_degradationEcoliDNA * ( ( t > 7.5*60 ) && ( t < 15*60 ) ); % molecules/cell/s
% % %     dNTPproduction = model.k_degradationEcoliDNA * ( ( t > 7.5*60 ) && ( t < 15*60 ) ); % molecules/cell/s
% %     % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % Host genome degradation
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % effective fraction of host genome remaining
    hostGenome = C(4 * model.n_gene + 9 + 1);
    
    if t < 7.5*60 
        
        hostGenomeDegrade = 0;
        
    elseif hostGenome > ( 0.15 * model.NhostGenomePerCell) ; % if there is sufficient host genome left
        
        hostGenomeDegrade = model.k_degradationEcoliDNA; % dNTP molecules/cell/s
        
    else
        
        % case when after 7.5 min but host genome is almost gone
        hostGenomeDegrade = 0; % fix this number?
        
    end
    
    % see below for dHostGenomedt
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % dNTP quantification
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    dntpProd = hostGenomeDegrade; % dNTP molecules/cell/s
    
%     % FBA Limits if virus draws more than is feasible
%     if viralLimits{1}(4 * model.n_gene + 8 + 1)
%         disp('limit to be evaled for dNTPs')
%     end
    
    
    if viralLimits{1}(4 * model.n_gene + 8 + 1) && (dntpProd > viralLimits{2}(4 * model.n_gene + 8 + 1))
        dntpProd = viralLimits{2}(4 * model.n_gene + 8 + 1);
%         
    end % limit viral flux and over limit
    hostGenomeDegrade = dntpProd;
    
    
    dHostGenomedt = - dntpProd / model.lengthHostGenome / 2;
    %
    %
    %  d(dNTP)                          2 * [bp of dNTP] * k_repT7RNA * [rep]
    % -------- = dNTPproduction  -  -------------------------------------- 
    %   dt                                  ([bp of dNTP] + K_mDNAelongat)
    %
    
    dntpCons = 2 * dnaProd/40 * model.lengthT7;
    
    ddNTPdt = dntpProd - dntpCons; % molecules/cell/s
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    
    
    
    
    
    
    % Rate Calculations for Procapsid
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Calculates rate of procapsid formation, as described
    % in Endy1997(Biotech&Bioeng). Procapsid formation begins to occur above
    % the threshold concentration of CnucleationT7procapsid using the
    % stoichiometry specified in Endy1997(Biotech&Bioeng) pg. 388.
    
    % Checks to see if threshold concentration of gp10A is reached
    if Cgp(model.genomeT7.gp10A.n_gene) > model.CnucleationT7procapsid
        % Once the threshold concentration of gp10A is reached then the rate of
        % procapsid formation, first term in equation 8 in Endy1997(Biotech&Bioeng)
        %
        %   k_assmblyT7procapsid * [ gp10A ^ nProcapsid ]
        %  ----------------------------------------------- = formationProCapsid
        %                    N_C
        %
        formationProCapsid = model.k_assemblyT7procapsid * Cgp(model.genomeT7.gp10A.n_gene) ^ (model.nProcapsid+1) / ...
            model.genomeT7.gp10A.stoichPhage;
    else
        % No procapsid formation occurse before the threshold concentration of gp10A is reached
        formationProCapsid = 0;
    end % procapsid nucleation
    
    % The net rate of procapsid concentration change is calculated using
    % equation 8 in Endy1997(Biotech&Bioeng)
    %
    %  d(PC)
    % ------ = formationProCapsid - ( k_packagingT7DNA ) * P_R
    %   dt
    %
    procapsidProd = formationProCapsid;
    
    % FBA Limits if virus draws more than is feasible
    if viralLimits{1}(3*model.n_gene+6) && (procapsidProd > viralLimits{2}(3*model.n_gene+6))
        procapsidProd = viralLimits{2}(3*model.n_gene+6);
    end % limit viral flux and over limit
    
%     procapsidCons = model.k_packagingT7DNA * CrateLimitingDNApackage;
    procapsidCons = phageProd / model.Nprogeney;
    
    dProCapsiddt = procapsidProd - procapsidCons;
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    
    

    
    
    
    % Combine the rates of change for function output
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    % placeholder (previously used for mRNA), gene products, T7 DNA,
    % T7 procapsid, T7 phage, dNTPs, and mRNAs for output and production
    % rates for FBA integration
    dCdt = [dGPdt; dDNAdt; dProCapsiddt; dPhagedt; ddNTPdt; dmRNAdt;...
        gpProd; dnaProd; procapsidProd; dPhagedt; dntpProd; mrnaProd; hostGenomeDegrade; dHostGenomedt];
    
    varargout = {dCdt};
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
else % identifiers genreation
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    [ identifiers ] = produceIdentifiers(model);

    varargout = {identifiers};
end

end % function ODEs_T7vEB1


% Subfunctions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% Calculate the translocation into host
function X = translocation(t,model)
    
    % translocation.m - Calculates the translocation of phage DNA (X in nt) 
    % at a given time (t in sec)
    %
    % Programmer:
    %   Elsa Birch - ewbirch@stanford.edu
    % 
    % Documentation:
    %   documentationT7vEB1.pdf (generated by documentationT7vEB1.tex)
    %
    % The code below calculates X(t) acording to equation (1) of the
    % ducumentation.
    %
    %
    % NOTE: possibly reorder comparisons for faster performance

    if t < model.tDelayInitial % initial delay
        X = 0;
    else 
        if t < model.tStartEcoliRNAP % basal translocation
            X = ( t - model.tDelayInitial ) * model.k_translocationT7(1); % nt
        else  
            if t < model.tStartT7RNAP % translocation mediated by E. coli RNA polymerase
                X = model.genomeT7.(model.switchTranslocation{1}).ntEnd + ...
                    ( t - model.tStartEcoliRNAP ) * model.k_translocationT7(2); % nt
            else
                if t < model.tFinishTranslocation % translocation mediated by T7 RNA polymerase
                    X = model.genomeT7.(model.switchTranslocation{2}).ntEnd + ...
                    ( t - model.tStartT7RNAP ) * model.k_translocationT7(3); % nt
                else % translocation finished
                    X = model.lengthT7;
                end % t < tFinishTranslocation  
            end % t < tStartT7RNAP
        end % t < tStartEcoliRNAP
    end % t < tDelayInitial
    X = floor(X);
end % translocation

% Calculate mRNA rates
function [dmRNAdt,capFlags,mRNAprod] = dmRNAdtCalculator(t,X,Cgp,CmRNA,model,viralLimits)

    % dmRNAdtCalculator- Calculates rate of change of mRNA concentration
    %
    % Programmers:
    %   Elsa Birch - ewbirch@stanford.edu
    %   Nicholas Ruggero - nruggero@stanford.edu
    %
    % Documentation:
    %   References to Endy1997B&B?, Endy1997Thesis
    %

    % Initalizations and counts to use later
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % Initialize rate array for dmRNA/dt
    mRNAprod = zeros(model.n_mRNA,1);
    capFlags = zeros(model.n_mRNA,2); % for debugging, flags to output when polymerase density capped at maximum 

    % rename count of promotors and 
    n_EcRNAP = model.n_rightEcRNAPpro;
    n_T7RNAP = model.n_rightT7RNAPpro;

    % determine number of each type in host cell at this time
    n_insertedEcRNAPpro = nnz( ((model.rightEcRNAPpro{3}(:,1)) < X)  ); % counts number of inserted EcRNAP promoters
    n_insertedT7RNAPpro = nnz( (model.rightT7RNAPpro{3} < X)  );

    % % determine which T7RNAP promotors are downstream of terminator and
    % % polymerase orrigionating from them are therefore not subject to termination
    % %
    % % Logical: 1 if promotor downstream of TT7 ntTerm
    % %          0 if promotor at or upstream of TT7 ntTerm
    % downstreamOfTermT7RNAPpro = (model.rightT7RNAPpro{3} > model.genomeT7.TT7.ntTerm );

    % all mRNA lengths, may not be used
    % % mRNALength = cell2mat(model.mRNA(:,4)) - cell2mat(model.mRNA(:,3)); 

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [Initalizations and counts to use later]


    % Determine active concentration of each ploymerase at this time t
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

    % EcRNAP
    %
    % polymerase = 0; % Sets the activePolymeraseConcentration function to calculate the free EcRNAP
    % false indicates that the calculation is performed for EcRNAP
    activeEcRNAP =  activePolymeraseConcentration(false,t,X,Cgp,model);

    % T7RNAP
    % polymerase = 1; % Sets the activePolymeraseConcentration function to calculate the free T7 RNAP
    % Only use fuction activePolymeraseConcentration.m which accesses
    % concentration of gp1, etc once this gene product is produced, otherwise
    % errors may result from function call.
    if Cgp(model.genomeT7.gp1.n_gene) > 0

        % Case when conentration of gp1 is nonzero
        % trueindicates that the calculation is performed for T7RNAP
        freeT7RNAP = activePolymeraseConcentration(true,t,X,Cgp,model);

    else
        % Case when conentration of gp1 is zero
        freeT7RNAP = 0; %

    end % if gp1 nonzero

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [Determine active concentration of each ploymerase at this time t]


    % (Prepare to calculate polymerase densities from each promotor at this time t)
    % First caclulate weighted length of transcript from each promotor, Le(i)
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % Initialize vectors of transcript length
    lengthTranscriptEcRNAP = zeros(n_EcRNAP, 1);
    lengthTranscriptT7RNAP = zeros(n_T7RNAP, 1);

    % Calculate Le for all inserted E. coli RNAP poromotors, k
    for k = 1:n_insertedEcRNAPpro
        % transcript length capped at TE for E. coli RNAP
        if X < model.genomeT7.TE.ntTerm

            % Before TE enters cell, assume transcript extends from last
            % nucleotide of promotor to last nucleotide in host cell
            %
            % Le = X(t) - EcRNAPpro_ntEnd
            lengthTranscriptEcRNAP(k) = X - model.rightEcRNAPpro{3}(k); % nt

        else
            % Once TE enters cell, assume transcript extends from last
            % nucleotide of promotor to nucleotide of TE termination
            %
            % Le = TE_ntTerm - EcRNAPpro_ntEnd
            lengthTranscriptEcRNAP(k) = model.genomeT7.TE.ntTerm - model.rightEcRNAPpro{3}(k); % nt

        end % if TE in cell
    end % for all inserted E. coli RNAP poromotors, k

    % Calculate Le for all inserted T7 RNAP poromotors, k
    for k = 1:n_insertedT7RNAPpro
        % transcript length capped at T7 for RNAP for promotors upstream of 
        % terminator, but with readthrough length accounted for.
        if X < model.genomeT7.TT7.ntTerm
            % Before TT7 enters cell, assume transcript extends from last
            % nucleotide of promotor to last nucleotide in host cell
            %
            % Le = X(t) - T7RNAPpro_ntEnd
            lengthTranscriptT7RNAP(k) = X - model.rightT7RNAPpro{3}(k); % nt

        elseif model.rightT7RNAPpro{2}(k) > model.genomeT7.TT7.ntTerm
            % TT7 has entered cell, but promotor k is downstream of terminator TT7, 
            % so assume transcript extends from last nucleotide of promotor to 
            % last nucleotide in host cell
            %
            % Le = X(t) - T7RNAPpro_ntEnd
            lengthTranscriptT7RNAP(k) = X - model.rightT7RNAPpro{3}(k); % nt

        else
            % TT7 has entered cell, and promotor k is upstream of terminator TT7, 
            % so assume transcript extends from last nucleotide of promotor to 
            % last nucleotide in host cell plus a contribution from readthrough
            % length weighted by the fraction of polymerase that continue
            % transcribing after TT7
            %
            % Le = (TT7_ntTerm - T7RNAPpro_ntEnd) + (1 - TT7_terminatorStrength)*(X(t) - TT7_ntTerm)
            lengthTranscriptT7RNAP(k) = model.genomeT7.TT7.ntTerm - model.rightT7RNAPpro{3}(k) ...
                                        + ( 1 - model.genomeT7.TT7.terminatorStrength) ...
                                        * ( X - model.genomeT7.TT7.ntTerm ); % nt
        end % if TT7 is in the cell
    end % for all inserted T7 RNAP poromotors, k

    % note that in lengthTranscriptEcRNAP, and lengthTranscriptT7RNAP the
    % lengths corresponding to promotors that are not yet in the host cell are
    % still zero.
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [First caclulate weighted length of transcript from each promotor, Le(i)]



    % Calculate polymerase densities from each promotor at this time t: Seb,
    % Sea, Stb, Sta
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % initialize arrays for densities of RNAP
    densityEcRNAPbefore = zeros(n_EcRNAP, 1);
    densityT7RNAPafter = zeros(n_T7RNAP, 1);
    densityT7RNAPbefore = zeros(n_T7RNAP, 1);
    % note that it is assumed here that no RNAP readthrough of TE, so Sea = 0
    % for all EcRNAP promotors

    % If any EcRNAP promotors inserted, calculate polymerase density from them
    % note that this calcualtion is occuring for each promotor, acomplished
    % using matrix elementwise operations so no for loop is needed
    if ~isempty(n_insertedEcRNAPpro)

        % First calcualte sum of protduct wieghted length and promotor strength for
        % all inserted EcRNAP promotors
        % __
        % \  Le_i * SA_i = sumInsertedEcRNAPproStrengthLength
        % /_ 
        %  i
        sumInsertedEcRNAPproStrengthLength = sum(lengthTranscriptEcRNAP(1:n_insertedEcRNAPpro) ...
                              .* (model.rightEcRNAPpro{4}(1:n_insertedEcRNAPpro,1))) + model.lengthHostGenome*0.25; % nt

        % Calculate polymerase density from each promotor
        % The density, Seb is calculated as
        %
        % Seb = [EcRNAP] * SA / sumInsertedEcRNAPproStrengthLength
        densityEcRNAPbefore(1:n_insertedEcRNAPpro) = activeEcRNAP ...
                        * (model.rightEcRNAPpro{4}(1:n_insertedEcRNAPpro,1)) ...
                        / sumInsertedEcRNAPproStrengthLength ; % should be RNAP molecule /nt

    end % if any EcRNAP promotors in cell


    % if any T7 promotors inserted, calculate polymerase density from them
    % Note that this calcualtion is occuring for each promotor, acomplished
    % using matrix elementwise operations so no for loop is needed
    if ~isempty(n_insertedT7RNAPpro)

        % First calcualte sum of protduct wieghted length and promotor strength for
        % all inserted T7RNAP promotors
        % __
        % \  Le_i * Sphi_i = sumInsertedT7RNAPproStrengthLength
        % /_ 
        %  i
        sumInsertedT7RNAPproStrengthLength = sum(lengthTranscriptT7RNAP(1:n_insertedT7RNAPpro) ...
                                  .* model.rightT7RNAPpro{4}(1:n_insertedT7RNAPpro,1) ); % nt

        % Calculate polymerase density from each promotor
        % The density, Stb is calculated as
        %
        % Stb = [T7RNAP] * Sphi / sumInsertedEcRNAPproStrengthLength 
        densityT7RNAPbefore(1:n_insertedT7RNAPpro) = freeT7RNAP ...
                            * model.rightT7RNAPpro{4}(1:n_insertedT7RNAPpro,1) ...
                            / sumInsertedT7RNAPproStrengthLength ; % should be RNAP molecule /nt

        % if terminator cacluate density afterward using the terminator
        % strength (fraction that terminate at TT7)
        % Note that the first term in the product is a logical array containing
        % 1 when the corresponding promotor (by index) is upstream of TT7, and
        % 0 when the promotor is after TT7 on the genome. This serves to assure
        % that the readthrough density calculated for promotors after the
        % terminator is 0.
        %
        % Sta = (1 - TT7_terminatorStrength) * [T7RNAP] * Sphi / sumInsertedEcRNAPproStrengthLength
        %     = (1 - TT7_terminatorStrength) * Stb
        densityT7RNAPafter = ((model.rightT7RNAPpro{3}(:,1)) < model.genomeT7.TT7.ntTerm) ...
                          .* (densityT7RNAPbefore) * ( 1 - model.genomeT7.TT7.terminatorStrength); % should be RNAP molecule /nt

    end % if any T7RNAP promotors in cell
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ Calculate polymerase densities from each promotor, Sea, Stb, Sta ]



    % (Prepare to calculate production rate for each mRNA)
    % Determine which promotors contribute to this mRNA, and sum RNAP denisties
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % for all mRNA (as defined in model.mRNA)
    for k = 1:model.n_mRNA
        % Transcription is only nonzero when full coding region (last nt) of mRNA is in host
        if model.mRNA{k,4} < X

            % EcRNAP contributions to transcription, calculate Se
            % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

    %         % Find E. coli promotors that are upstream of mRNA k start 
    %         %
    %         % Logical: 1 if promotor upstream of mRNA start 
    %         %          0 if mRNA starts before promotor
    %         mRNAdownstreamEcRNAPpro = (model.rightEcRNAPpro{2}(:,1)) <= model.mRNA{k,3};

    %         % Check if mRNA k exceeds terminator, therefore promotors contribute 
    %         % only readthrough polymerase density 
    %         %
    %         % Logical: 0 if mRNA_ntEnd < TE_ntTerm  (mRNA ends downstream of TE)
    %         %          1 if mRNA_ntEnd >= TE_ntTerm (mRNA ends at or upstream of TE)
    %         mRNAupstreamEcRNAPterm =  model.genomeT7.TE.ntTerm >= model.mRNA{k,4};

            % Sum E. coli RNAP density contribution of promotors upstream of mRNA k
            if model.mRNAupstreamEcRNAPterm(k)

                % Total E. coli RNAP density over coding region of mRNA k 
                % Note: elementwise multiplication by logical mRNAdownstreamEcRNAPpro
                % enforces condition that promotor is upstream of mRNA k
                %      __
                % Se = \  Seb_i  , for i over promotors upstream of this mRNA k
                %      /_                as enforced by 0 logical value of mRNAdownstreamEcRNAPpro  
                %       i                for promotors after mRNA start multiplied elementwise
                %
                %  ( note Seb_i 0 for any uninserted promotors, see above )
                Se = sum(model.mRNAdownstreamEcRNAPpro{k} .* densityEcRNAPbefore);

            else
                % Case when mRNA ends after TE
                % Assuming 100% termination at TE, no readthrough polymerase
                % density  over coding region of mRNA k
                Se = 0;

            end % if mRNA is upstream of TE

            % Note that in WT T7 genome there are no E coli RNAP promotors
            % after the terminator TE, this the two catches above are sufficent

            % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::[ EcRNAP contributions to transcription, calculate Se ]


            % T7 RNAP contributions to transcription, calulcate St
            % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

    %         % Find T7 promotors that are upstream of mRNA k start 
    %         %
    %         % Logical: 1 if promotor upstream of mRNA k start 
    %         %          0 if mRNA starts before promotor
    %         mRNAdownstreamT7RNAPpro = (model.rightT7RNAPpro{2}(:,1)) <= model.mRNA{k,3};

    %         % Check if mRNA k exists entirely before terminator 
    %         %
    %         % Logical: 0 if mRNA_ntEnd < TT7_ntTerm  (mRNA ends downstream of TT7)
    %         %          1 if mRNA_ntEnd >= TT7_ntTerm (mRNA ends at or upstream of TT7)
    %         mRNAupstreamT7RNAPterm =  model.genomeT7.TT7.ntTerm >= model.mRNA{k,4};
    %         
    %         % Check if mRNA k exists entirely after terminator
    %         %
    %         % Logical: 0 if mRNA_ntStart < TT7_ntTerm  (mRNA starts upstream of TT7)
    %         %          1 if mRNA_ntStart > TT7_ntTerm (mRNA starts downstream of TT7)
    %         mRNAdownstreamT7RNAPterm =  model.genomeT7.TT7.ntTerm < model.mRNA{k,3};

            % Sum T7 RNAP density contribution of promotors upstream of mRNA k,
            % (mothod choice based on position of mRNA with respect to terminator)
            if model.mRNAupstreamT7RNAPterm(k)

                % This is the case when the entire mRNA is before the
                % terminator TT7, so only before termination RNAP densities are
                % relevent to transcription
                %      __
                % St = \  Stb_i  , for i over promotors upstream of this mRNA k
                %      /_                as enforced by 0 logical value of mRNAdownstreamT7RNAPpro  
                %       i                for promotors after mRNA start multiplied elementwise
                %
                %  ( note Stb_i 0 for any uninserted promotors, see above )
                St = sum( model.mRNAdownstreamT7RNAPpro{k} .* densityT7RNAPbefore );

            elseif model.mRNAdownstreamT7RNAPterm(k)

                % This is the case when the entire mRNA is after the
                % terminator TT7, so some readtrhough (after) RNAP densities are
                % relevent to transcription, but also some densities of RNAP
                % orrigionating from promotors downstream of TT7 (before)
                % contribute
                %      __
                % St = \  Stb_i  , for i over promotors upstream of this mRNA k
                %      /_                as enforced by 0 logical value of mRNAdownstreamT7RNAPpro  
                %       i                for promotors after mRNA start multiplied elementwise
                %                        and (promotors) downstream of term TT7
                %                        as enforced by 0 logical value of downstreamOfTermT7RNAPpro
                %       __               for promotors upstream of TT7 multiplied elementwise.
                %     + \  Sta_j  , for j over promotors upstream of this mRNA k
                %       /_                as enforced by 0 logical value of mRNAdownstreamT7RNAPpro  
                %        j                for promotors after mRNA start multiplied elementwise
                %                         and (promotors) upstream of term TT7
                %                         as enforced by 0 logical value of not(downstreamOfTermT7RNAPpro)
                %                         for promotors downstream of TT7 multiplied elementwise.
                %
                %  ( note Stb_i and Sta_j 0 for any uninserted promotors, see above )
                St = sum( model.downstreamOfTermT7RNAPpro .* model.mRNAdownstreamT7RNAPpro{k} .* densityT7RNAPbefore ) ...
                     + sum( ~model.downstreamOfTermT7RNAPpro .* model.mRNAdownstreamT7RNAPpro{k} .* densityT7RNAPafter );

            else
                % This is the case when the entire mRNA spanse the
                % terminator TT7, so only readthrough (after) RNAP densities are
                % relevent to transcription
                %      __
                % St = \  Sta_i  , for i over promotors upstream of this mRNA k
                %      /_                as enforced by 0 logical value of mRNAdownstreamT7RNAPpro  
                %       i                for promotors after mRNA start multiplied elementwise
                %
                %  ( note Sta_i 0 for any uninserted promotors, see above )
                St = sum( mRNAdownstreamT7RNAPpro .* densityT7RNAPafter );

            end % mRNAupstreamT7RNAPterm
            % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::[ T7 RNAP contributions to transcription Calculate St ]



            % Cap polymerase densities
            % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

            % Set E coli RNA polymerase density at maximum value, if exceeded
            if Se > model.spacingRNAploymerase
                Se = model.spacingRNAploymerase;
                % capFlags(k,2) = 1;
            end % is Ec RNAP density exceeds maximum

            % Set T7 RNA polymerase density at maximum value, if exceeded
            if St > model.spacingRNAploymerase
                St = model.spacingRNAploymerase;
                % capFlags(k,1) = 1;
            end % is T7 RNAP density exceeds maximum 

            % :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::[ Cap polymerase densities ]

            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -[Determine which promotors contribute to this mRNA, and sum RNAP denisties]



            % mRNA ODE evaluation
            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            % Evaluate mRNA ODE as in Endy1997Thesis eqn (x), for each mRNA k
            %
            % d(mRNA)
            % ------ = (k_transcriptionEcoli * Se) + (k_transcriptionT7 * St) - (k_decayT7mRNA * [mRNA])
            %   dt   
            %       

            mRNAprod(k) = model.k_transcriptionEcoli*Se + model.k_transcriptionT7*St;

    %         % FBA Limits if virus draws more than is feasible
    %         if viralLimits{1}( 3*model.n_gene+8 + k) && (mRNAprod(k) > viralLimits{2}( 3*model.n_gene + 8 + k))
    %             mRNAprod(k) = viralLimits{2}( 3 * model.n_gene + 8 + k);
    %         end % limit viral flux and over limit


    %         mRNAcons(k) = CmRNA(k) * model.k_decayT7mRNA;



            % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        end % if mRNA start nt is in cell

    end % for all mRNA k

    % apply limits from FBA as needed
    mRNAprodLimitFlags = viralLimits{1}( 3*model.n_gene+8 + 1 : 4*model.n_gene+8 );
    if any( mRNAprodLimitFlags )

        mRNAprodLimits = viralLimits{2}( 3*model.n_gene+8 + 1 : 4*model.n_gene+8 );
    %     size( mRNAprod ) 
    %     size( mRNAprodLimits )
        logicalMRNAoverLimit = mRNAprodLimitFlags & ( mRNAprod > mRNAprodLimits );
        mRNAprod( logicalMRNAoverLimit ) = mRNAprodLimits( logicalMRNAoverLimit );

    end


    % mRNA ODE evaluation
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    % Evaluate mRNA ODE as in Endy1997Thesis eqn (x), for each mRNA k
    %
    % d(mRNA)
    % ------ = (k_transcriptionEcoli * Se) + (k_transcriptionT7 * St) - (k_decayT7mRNA * [mRNA])
    %   dt   
    %

    mRNAcons = CmRNA * model.k_decayT7mRNA;

    dmRNAdt = mRNAprod - mRNAcons;% should be molecules/cell

end % dmRNAdtCalculator

% Calculate Gene Product rates
function [dGPdt,gpProd] = dGPdtCalculator(Cgp,CmRNA,~,X,model,viralLimits, ratePhageAssembly)

    % dGPdtCalculator.m - Calculates rate of change of gene product
    % concentrations
    %
    % Programmers:
    %   Elsa Birch - ewbirch@stanford.edu
    %   Nicholas Ruggero - nruggero@stanford.edu
    %
    % Documentation:
    %   References to Endy1997B&B?, Endy1997Thesis
    %
    % Dependancies:
    %   Called by: ODEs_t7vEB1.m
    %   Calls: proteinsToAssembly.m
    %
    % Caclulates three things related to the rate of change
    % in the concentrations of gene products from the T7 genome.
    %
    % Rate of change of concentrations of gene products from the T7 genome is calculated.

    % Initialize rate array for dgp/dt
    dGPdt = zeros(model.n_gene,1);
    gpProd = zeros(model.n_gene,1);
    gpCons = zeros(model.n_gene,1);

    % ribosome avalibility calculations
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    % calculate total nts of mRNA
    ntTotalmRNA = sum( (CmRNA .* (model.mRNAcompact{4} - model.mRNAcompact{3})) ); % nt

    if ntTotalmRNA > 0

        densityRibosomes = model.n_ribosome / ntTotalmRNA ; % ribo/nt

        if densityRibosomes > (1/model.spacingRibosome)

            densityRibosomes = (1/model.spacingRibosome);

        end

    else
        densityRibosomes = 0;

    end

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

    % [cRateLimitingDNApackage, rateLimitPhageAssembly ] = determineT7assemblyRateLimit(C, model);

    % For all phage T7 gene products gp k, prepare parameters for and complete
    % cauculation of gp ODE
    for k = 1 : model.n_gene



        % Only allow translation once coding region for protien has
        % entered, this allows gene product production before whole coding
        % region for gene has been translocated into cell. Justified by
        % Endy1997Thesis code
        if (model.gene{2}(k,1) < X)

            % Find total cellular concentration of mRNAs containing gene
            mRNAcontainingGP = model.gene{6}(k,1);
            %
            % [mRNA of gp] = 
            totalGPmRNA = sum(CmRNA(mRNAcontainingGP)); % mol/cell

            % Set concentration of gp3.5 mRNA equal to that of gp2.8 mRNA as in
            % Endy1997 Thesis
            % % if k == model.genomeT7.gp3o5.n_gene
            % %    mRNAcontainingGP = model.gene{k-2,6};
            % %    totalGPmRNA = sum(CmRNA(mRNAcontainingGP));
            % % end % this gene is 3.5

            % Find length of gene in nt, for use later
            %
            % geneLength = gp_ntEnd - gp_ntStart
            % % geneLength = model.gene{k,3} - model.gene{k,2};

            % Call proteinsToAssembly function to calcualte the number of proteins
            % which are used in assembly of the T7 procapsid, T_gp which is output
            % here as an array where T_gp is negative when gp is being consumed by
            % phage particle assembly
            % 
            % T_gp = gpToAssembly
            if model.gene{4}(k) > 0
                gpToAssembly = proteinsToAssembly(k, Cgp, model, ratePhageAssembly);
            else
                gpToAssembly = 0 ;
            end


            % Evaluate gp rate ODE from Endy1997 Thesis
            % equation (x), for each gp, k
            %
            %  d(gp)    k_translation * [mRNA of gp]
            % ------ = -----------------------------  -  k_decayProtien * [gp] + T_gp 
            %   dt              S_ribosome
            %




            gpProd(k) = totalGPmRNA * model.k_translation * densityRibosomes;
            gpCons(k) = - model.k_decayProtien * Cgp(k) + gpToAssembly;

            % FBA Limits if virus draws more than is feasible
            if viralLimits{1}( 2 * model.n_gene + 4 + k) && (gpProd(k) > viralLimits{2}( 2 * model.n_gene + 4 + k))
                gpProd(k) = viralLimits{2}( 2 * model.n_gene + 4 + k);
            end % limit viral flux and over limit


            dGPdt(k) = gpProd(k) + gpCons(k); % molecules/cell



        end % if gene in cell
    end % for k over all gp

end % dGPdtCalculator

% Calculate binding equilibrium
function [Cbound] = bindingEquilib(Ctotals,Keq,vol_Ecoli,t)


% ODEs_t7vEB1 - ODEs for implementation of t7v1 from Endy1997(Biotech&Bioeng)
%
% Programmers:
%   Elsa Birch - ewbirch@stanford.edu
%   Nicholas Ruggero - nruggero@stanford.edu
%
% Documentation:
%   References to Endy1997B&B, Endy1997Thesis
%
% Dependancies:
%   Called by: setup_t7vEB1.m
%   Calls: - translocation.m
%          - dmRNAdtCalculator.m (calls, activePolymeraseConcentration.m (calls bindingEquilib.m))
%          - dGPdtCalculator.m (calls proteinsToAssembly.m)
%
% Chemical kinetic (ordinary differential equation model) of bacteriophage
% T7 replication in host E. coli
%
% Notes:
%   - Assumes temperature is 30C
%   - Time unit used is 1 second
%   - T7 genetic elements named with decimal converted to 'o' in variable
%   names, eg 1.1A becomes 1o1A

    % C threshold
    Cthresh = 1;

    if any( Ctotals < Cthresh )
        
        % Case where concentration of either polymerase or inhibitor is
        % zero, therefore no binding occurs
        Cbound = 0;
        
    else
        
        % First convert concentrations
        N_A = 6.022e23; % Avagadro's number for unit conversion
        Ctotals = Ctotals / N_A / vol_Ecoli; % unit conversion from molecules/cell to M
   
        % Solve quadratic resulting from binding equilibrium equation and
        % mass balances using the quadratic formula as in Endy1997Thesis
        % code, using p1 & p2 as the protiens binding
        %
        %         1 (
        % bound = - ( ( [p1_total] + [p2_total] + 1/K_eq ) - 
        %         2 (
        %                  __________________________________________
        %            -    / ( [p1_total] + [p2_total] + 1/K_eq )^2 
        %               \/   
        %                          __________________________________   )
        %                           - 4 * [p1_total] * [p2_total] ]     )
        %                                                               )
        a = 1;
        b = (Ctotals(1) + Ctotals(2) + 1/Keq);
        c = Ctotals(1) * Ctotals(2);
        Cbound = (b - sqrt(b^2 - 4*a*c))/(2*a);
        
        
        
        % Warn if mass balance violated, since that is theoretically
        % possible using above equation
        if any(Cbound > Ctotals)
            disp('mass bal violated')
            Cbound 
            Ctotals
        end
        
        
        % This is the orrigional method implemented to solve for roots of
        % equilibrium binding quadratic, using logic to assure physicality
        % of bound concentration found
        %::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
        % % % Calculates roots of equation specified in the notes
        % % N_A = 6.022e23; % Avagadro's number for unit conversion
        % % Ctotals = Ctotals / N_A / vol_Ecoli; % unit conversion from molecules/cell to M
        % % 
        % % % Set up quadratic coefficients defining equilibrium
        % % A = Keq * Ctotals(1) * Ctotals(2); % zero order term
        % % B = - ( ( Ctotals(1) + Ctotals(2) ) * Keq + 1); % first order term
        % % eqQuadPolynomial = [ Keq B A ];
        % % eqQuadPolyRoots = roots( eqQuadPolynomial ); % roots
        % % 
        % % % Checks roots for existance & physical reasonability
        % % if isempty( eqQuadPolyRoots )
        % %     % if no roots found generate error and assume none bound
        % %     Cbound = 0;
        % %     disp(['no roots found for binding with Keq (', num2str(Keq),...
        % %           '), Ctotals (', num2str(Ctotals(1)), ', ',...
        % %           num2str(Ctotals(2)), ') at time t = ', num2str(t)]);
        % % else
        % %     % calculate number of physically feasible roots
        % %     n_feasible = sum( ( eqQuadPolyRoots > 0 ) & ( eqQuadPolyRoots < min(Ctotals) ));
        % % 
        % %     switch n_feasible
        % %         case 0
        % %             % if no physically feasible roots found generate error and assume none bound
        % %             Cbound = 0;
        % %             disp(['no physically feasible found for binding with Keq (', ...
        % %                    num2str(Keq),'), Ctotals (', num2str(Ctotals(1)), ', ',...
        % %                    num2str(Ctotals(2)), ') at time t = ', num2str(t)]);
        % %         case 1
        % %             % use one physically feasible root as concentration bound
        % %             Cbound = eqQuadPolyRoots( find( ( eqQuadPolyRoots > 0 ) & ( eqQuadPolyRoots < min(Ctotals) ) ) );
        % %         case 2
        % %             % if 2 physically feasible roots found generate error and assume max bound
        % %             % NOTE: Junk checks for one root being essentially zero
        % %             [Cbound, junk] = max(eqQuadPolyRoots);
        % %             if junk > sqrt(eps)                    
        % %                 disp(['2 physically feasible found for binding with Keq (', ...
        % %                        num2str(Keq),'), Ctotals (', num2str(Ctotals(1)), ', ',...
        % %                        num2str(Ctotals(2)), ') at time t = ', num2str(t), ...
        % %                        ', maximum used of ', num2str(eqQuadPolyRoots(1)),', '...
        % %                        num2str(eqQuadPolyRoots(2))]);
        % %             end % end of check if junk is > eps
        % %     end
        % % end
        %::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
        
        % Convert units back for use in model
        Cbound = Cbound * N_A * vol_Ecoli; % Convert back to molecules/cell
        
       
    end % if both Ctotals nonzero
end % bindingEquilib

% Calculate phage componenet that limits phage assembly
function [cRateLimitingDNApackage, rateLimitPhageAssembly ] = determineT7assemblyRateLimit(c, model)
    % Determine the rate limiting species for phage formation

    % Handle DNA and procapsids
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % There was no equation given but the limiting
    % species was listed as either DNA or procapsids so assumed to be
    % lesser of these at time t
    %
    % P_R = minimum{[PC], [T7DNA]}
    %

    cRateLimitingDNApackage = min(c((model.n_gene + 1):(model.n_gene + 2))); % molecules/cell
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ Handle DNA and procapsids ]



    % Handle all other GP
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    % GP concentrations
    cGP = c(1:model.n_gene); % molecules/cell

    % Stoichiometric coeffs in phage particles for all GP
    stoichGP = (model.gene{4}(:,1));

    % Indicies of GPs included in phage particle
    indGPinPhage = find(stoichGP);

    % Devide all GP concentrations by stoich phage coeff to determine phage
    % possible from individual gene products at that time
    %
    %                           [gp_i]
    % phageAssemblyLimit_i = -------------- = [phage_i]
    %                         gp_i / phage
    %
    % accomplished using elementwise division
    gpForLimit = cGP(indGPinPhage) ./ stoichGP(indGPinPhage); % phage molecules/cell

    % minimum phage assembly rate limited by GPs
    %     cMinGP = min(phageAssemblyLimit); % molecules/cell

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ Handle all other GP ]

    % Compare two potential limiters of phage assembly rate
    rateLimitPhageAssembly = min([gpForLimit; c((model.n_gene + 1):(model.n_gene + 2))]);
    
end % determineT7assemblyRateLimit

% Calculate the amount of each protien going towards assembled phage
function [ gpToAssembly] = proteinsToAssembly(k, Cgp, model, ratePhageAssembly)

% proteinsToAssembly.m - Calculates the amount of gp9 and gp10A which
% are used in the T7 procapsid assembly
%
% Documentation:
%   References to Endy1997B&B?, Endy1997Thesis
%
% Dependancies:
%   Called by: dGPdtCalculator.m
%   Calls: 
%
% Programmers:
%   Elsa Birch - ewbirch@stanford.edu

    % Calculates the gene product used in procapsid/phage assembly, eqn (9a,b,c)
    % and determines the limiting species for T7 assembly
    
    % first determine which GP, or PC or DNA limits the rate of phage
    % assembly currently
%     CrateLimitingDNApackage = determineT7assemblyRateLimit(C, model);
    
    % Calculates the amount of gene products that are being incorperated in
    % procapsids and viral progeney particles, T_gp
    switch model.gene{1}{k}
        case 'gp9'
            if Cgp(model.genomeT7.gp10A.n_gene) > model.CnucleationT7procapsid
                % In the case that the threshhold concentration of major
                % capsid protien 10A, C_N,
                % Gene product 9 scaffold protiens used (not reused) in
                % procapsid assembly calculated as:
                %
                %                  k_assmblyT7procapsid * [ gp10A ^ nProcapsid ]
                % T_gp9 = -(N_s)( ----------------------------------------- )
                %                                 N_C
                %
                % Here, N_s is (scaffolds gp9 / procapsid) and 
                %       N_C is (major capsid protien gp10A / procapsid)
                gpToAssembly = - model.genomeT7.gp9.stoichProcapsid * ...
                               model.k_assemblyT7procapsid * ...
                               Cgp(model.genomeT7.gp10A.n_gene) ^ (model.nProcapsid+1) / ...
                               model.genomeT7.gp10A.stoichPhage;
            else
                % Concentration gp10A below critical procapsid nucleation,
                % therefore no gp9 used in procapsid assembly at this t
                %                   
                % T_gp9 = 0
                %
                gpToAssembly = 0;
            end % if procapsid nucleation concentration reached
                       
        case 'gp10A'
            
            if Cgp(model.genomeT7.gp10A.n_gene) > model.CnucleationT7procapsid
                % In the case that the threshhold concentration of major
                % capsid protien 10A, C_N, Gene product consumed in 
                % procapsid assembly calculated as:
                %                   
                % T_gp10 = - k_assmblyT7procapsid * [ gp10A ^ nProcapsid ]
                %
                gpToAssembly = - model.k_assemblyT7procapsid * ...
                               Cgp(model.genomeT7.gp10A.n_gene) ^ (model.nProcapsid+1);
            else
                % Concentration gp10A below critical procapsid nucleation,
                % therefore no gp10A used in procapsid assembly at this t
                %                   
                % T_gp10 = 0
                %
                gpToAssembly = 0;
            end % ifprocapsid nucleation concentration reached
            
        otherwise
            % Case for all other gene products
            % evaluate the general expression for use of phage protiens in
            % procapsid assembly, Endy1997Thesis eqn (3.5c)
            %
            % T_gpk = - (N_k) * ( N_progeney ) * ( k_packagingT7DNA ) * P_R
            %
            % Here N_k is the stoichiometric coefficient of gp k protiens
            %          to progeney assembled, 0 if not included in progeny
            %      N_progeney is efficieny of DNA or Procapsids (P_R)
            %                 resulting in functional progeney, from ENdy1997B&B
            
%             ratePhageAssembly
            
%             gpToAssembly = - model.gene{4}(k) * model.Nprogeney *...
%                            model.k_packagingT7DNA * CrateLimitingDNApackage;
                       
            gpToAssembly = - model.gene{4}(k) * ratePhageAssembly;
                       
                       
    end % switch for gp9, 10A or all others
    
end % proteinsToAssembly

% Calculate remaining active polymerase
function [activePolymerase] = activePolymeraseConcentration(polymerase,t,X,Cgp,model)

% activeEcRNAPConcentration.mat - Calculates the concentration of
% EcRNAP that is not inhibited by gp0.7 or gp2.0 and T7 RNAP that is not
% inhibited by gp3.5
%
% Programmers:
%   Elsa Birch - ewbirch@stanford.edu
%   Nicholas Ruggero - nruggero@stanford.edu
%
% Documentation:
%   References to Endy1997B&B, Endy1997Thesis
%
% Dependancies:
%   Called by: dmRNAdtCalculator.m 
%   Calls: - bindingEquilib.m

% If polymerase = 0 then the calculation is performed for EcRNAP
% If polymerase = 1 then the calculation is performed for T7 RNAP


% EcRNAP active polymerase determination
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
if polymerase == 0
    
    % If the end of Class I has been traslocated then the length of
    % Class I DNA inserted is the last nucleotide of Class I DNA. Otherwise
    % it is whatever is the last translocated nucleotide
    % % ntEndClass1 = model.genomeT7.TE.ntTerm;
    % % if X < ntEndClass1 
    % %     % % lengthInsertedClass1 = X;
    % %     lengthInsertedClass1 = ntEndClass1;
    % % else
    % %     lengthInsertedClass1 = ntEndClass1;
    % % end % X < ntEndClass1
    
    % First order EcRNAP decay with the calculated result in molecules/cell
    % Endy1997(Biotech&Bioeng) described in methods
    %
    % EcRNAP_total = EcRNAP_initial * e ^ ( - k_decayProtien * t )
    %
    totalEcRNAP = model.n_polymeraseEcoli * exp((-model.k_decayProtien)*t);
    
    % Note: 30% attributed to gp0o7 and 70% attributed to gp2o0 as stated
    % in Endy1997(Biotech&Bioeng)
    %
    gp0o7boundEcRNAP = bindingEquilib([Cgp(model.genomeT7.gp0o7.n_gene), (totalEcRNAP)], model.K_eqBindGp0o7, model.vol_Ecoli,t);
    gp2boundEcRNAP = bindingEquilib([Cgp(model.genomeT7.gp2.n_gene), ( totalEcRNAP)], model.K_eqBindGp2, model.vol_Ecoli,t);
    
    % Free EcRNAP as determined by being bound to niehter gp0o7 or gp2 , 
    % in Endy1997(Biotech&Bioeng) the caculation implied in methods is
    % subtracting both from total after a 30/70 split of possible
    % inhibition by either protien. However Endy1997Thesis code shows the
    % calculation of inhibition by both using the full total of EcRNAP then
    % the inhibitiion used as being the maximum of inhibitors
    %
    freeEcRNAP = (totalEcRNAP - max([gp2boundEcRNAP,gp0o7boundEcRNAP]));
%     freeEcRNAP = (totalEcRNAP - gp2boundEcRNAP - gp0o7boundEcRNAP);% molecules/cell
    
    % Used in past variation of inhibition calculation
    % % freeEcRNAP < 0 
    % %     freeEcRNAP = 0;
    % % end
    
    % Method inmplied in methods in Endy1997(Biotech&Bioeng)
    % % fracActiveEcRNAP = freeEcRNAP / totalEcRNAP; % Called percent fraction in text, see notes 
	% % activePolymerase = lengthInsertedClass1 * fracActiveEcRNAP / model.spacingRNAploymerase; % P_h in molecules/cell
    
    % Method inmplied in methods/code in Endy1997Thesis
    activePolymerase =  freeEcRNAP;
    
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ EcRNAP active polymerase determination ]
    

% T7RNAP active polymerase determination
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
elseif polymerase == 1
    % length of class II & III T7 DNA in the host
    % calculate 
    % % ntEndClass1 = model.genomeT7.phi1o1A.ntStart;
    % % lengthInsertedClass23 = X - ntEndClass1;
    
    % Binding of gp1 by gp 3.5 
    gp3o5boundT7RNAP = bindingEquilib( [Cgp(model.genomeT7.gp1.n_gene),Cgp(model.genomeT7.gp3o5.n_gene)],...
                        model.K_eqBindGp3o5, model.vol_Ecoli ,t); % Assumed from text, molecules/cell
    
    % Method inmplied in methods in Endy1997(Biotech&Bioeng)
    % % freeT7RNAP = (Cgp(model.genomeT7.gp1.n_gene) - gp3o5boundT7RNAP) / Cgp(model.genomeT7.gp1.n_gene);
	% % activePolymerase = lengthInsertedClass23 * freeT7RNAP / model.spacingRNAploymerase;% Assumed from text, molecules/cell
    
    % Method inmplied in code in Endy1997Thesis
    activePolymerase = (Cgp(model.genomeT7.gp1.n_gene) - gp3o5boundT7RNAP);
    

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [ T7RNAP active polymerase determination]
    
end % which polymerase activity is being calculated

end

% Produce cell array of identifying strings from information in model structure
function [ identifiers ] = produceIdentifiers(model)
    identifiers = repmat({''},(2*(2*model.n_gene+4)+2),2);
%     size(identifiers)
    for iMetRxn = 1:2*(2*model.n_gene+4)+1
        if iMetRxn <= model.n_gene% viral gp accumulation
            identifiers{iMetRxn,1} = model.gene{1}{iMetRxn,1};
            
        elseif model.n_gene+1 == iMetRxn % viral dna accumulation
            identifiers{iMetRxn,1} = 'genome';
            
        elseif model.n_gene+2 == iMetRxn % viral procapsid accumulation
            identifiers{iMetRxn,1} = 'procapsid';
            
        elseif model.n_gene+3 == iMetRxn % viral phage accumulation
            identifiers{iMetRxn,1} = 'T7';
            
        elseif model.n_gene+4 == iMetRxn % viral dntp accumulation
            identifiers{iMetRxn,1} = '';
            
        elseif model.n_gene+4+1 <= iMetRxn && iMetRxn <= 2*model.n_gene +4 % viral mrna accumulation
            % identifiers{iMetRxn,1} = ['mRNA', model.gene{iMetRxn - (model.n_gene+4)}];
            geneNameString = model.gene{1}{iMetRxn - (model.n_gene+4),1};
            identifiers{iMetRxn,1} = ['mRNA', geneNameString(3:end)];
            
        elseif 2*model.n_gene+4 + 1 <= iMetRxn && iMetRxn <= 2*model.n_gene+4 + model.n_gene % viral gp production
            identifiers{iMetRxn,2} = ['R_protein_',model.gene{1}{iMetRxn - (2*model.n_gene+4),1}];
            
        elseif 2*model.n_gene+4 + model.n_gene + 1 == iMetRxn % viral dna production
            identifiers{iMetRxn,2} = 'R_genome';
            
        elseif 2*model.n_gene+4 + model.n_gene + 2 == iMetRxn % viral procapsid production
            identifiers{iMetRxn,2} = 'R_procapsid';
            
        elseif 2*model.n_gene+4 + model.n_gene + 3 == iMetRxn % viral phage production
            identifiers{iMetRxn,2} = 'R_phage';
            
        elseif 2*model.n_gene+4 + model.n_gene + 4 == iMetRxn % viral dntp production
            identifiers{iMetRxn,2} = '';
            
        elseif 2*model.n_gene+4 + model.n_gene+4+1 <= iMetRxn && iMetRxn <= 2*model.n_gene+4 + 2*model.n_gene+4 % viral mrna production
            identifiers{iMetRxn,2} = ['R_mrna_',model.gene{1}{iMetRxn - (2*model.n_gene+4) - (model.n_gene+4),1}];
        elseif 4 * model.n_gene + 8 + 1 == iMetRxn
            identifiers{iMetRxn,2} = 'R_ECGENDEG';
        elseif 4 * model.n_gene + 9 + 1 == iMetRxn
            identifiers{iMetRxn,2} = ''; % this is the host genome
        end
    end
    
end


