% modelFromScratch.m - All constant and calculated from constant parameters
% for model are initalized here and stored in a structure called 'model'
%
% Programmers:
%   Elsa Birch - ewbirch@stanford.edu

function model = T7_37Growth_EB_VirusModelFromScratch(mu) %, k_assemblyT7procapsid, k_packagingT7DNA, ...
    % k_replicationT7DNA, fTrans, fRibo, fPoly, fTranscr, fSpace, nProcapsid)



% Constant Parameters
% SOME OF THESE ARE RECALCULATED LATER AS HOST MU DEPENDANT
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
model = struct(...
     ...% simulation parameters [Table I in Endy1997(Biotech&Bioeng)]
     ...% kinetic
        'k_translocationT7',        [ 141 , 55 , 400 ],... % bp/sec
        'k_transcriptionEcoli',     55,... % nt/sec/RNAP
        'k_transcriptionT7',        400,... % nt/sec/RNAP
        'k_decayT7mRNA',            0.001,...0 % /sec
        'k_translation',            42,... nt/sec/ribo 14,... % AA/sec/ribosome
        'k_decayProtien',           3.92e-5,... % /sec
        'k_degradationEcoliDNA',    32357,... % nucleotides/sec
        'k_replicationT7DNA',       475,... k_replicationT7DNA,...370,... % bp/sec/polymerase
        'k_packagingT7DNA',         0.983,...k_packagingT7DNA,... % /min
        'k_assemblyT7procapsid',    6.44e-16,...k_assemblyT7procapsid,...6.44e-16,... % /(number/cell)^3.78/min
     ...% binding
        'K_eqBindGp2',              5.0e7, ... %5.0e7,... % /M 
        'K_eqBindGp0o7',            5.5e6,... %5.5e6,... % /M
        'K_eqBindGp3o5',            1.087e7,... % /M
     ...% other
        'eta_terminationT7',        0.66,... % efficeincy
        'n_polymeraseEcoli',        1800,... % molecules/cell
        'K_mDNAelongation',      	8668,... % nucleotides/cell
        'CnucleationT7procapsid',   1000,...,... % molecules/cell
        'vol_Ecoli',                8e-16,... % L
        'Nprogeney',                1,... % T7 progeney per procapsid/DNA
        'NhostGenomePerCell',       1.84,... % number genome equaivalents
        'fractHostGenomeDegrade',   0.85,... % fraction degraded to be avalible to phage
        'lengthHostGenome',         4.655e6,...% bp
        'tTotalHostGenomeDegrade',  450, ...  % sec
        'nProcapsid',               4.78... nProcapsid...% 
        ); 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
if isempty(mu)
    % Growth rate
    model.hostGrowth = 1; % /h    ref endy thesis code
else
    
    model.hostGrowth = mu; % /h 
    
end
  
if model.hostGrowth < 0.6
    model.hostGrowth = 0.6;
end

fTrans = 1;
fRibo = 1;

% Calculated parameters
% Growth rate dependant simulation parameters 
% [Table I in Endy1997(Biotech&Bioeng)]
model.n_polymeraseEcoli = -103.59+model.hostGrowth*(548.977+model.hostGrowth*(-155.18+192.216*model.hostGrowth));
model.spacingRNAploymerase = 233/model.hostGrowth^2 + 73; % bp
% model.spacingRibosome = 82.5/model.hostGrowth + 145; % nt/ribosome
model.spacingRibosome = ((((-31.7460*model.hostGrowth+215.5556)*model.hostGrowth-519.4921)*model.hostGrowth+492.7778)*model.hostGrowth-72.0952 );%* fSpace; % nt/ribosome
% % model.spacingRibosome = model.spacingRibosome / 3; % aa/ribosome
c = 42.4 + 132.7*exp(-2.812 * model.hostGrowth); % min
d = 22.3 + 15.4 * exp(-1.16 * model.hostGrowth); % min
model.NhostGenomePerCell = 60 / model.hostGrowth / (c * log(2) ) * ( 2 ^ (( c + d ) * model.hostGrowth / 60 ) - 2 ^ ( model.hostGrowth / 60 ) );
model.k_translation = (102.9 - 55.5 * (model.hostGrowth) ^ (-0.3906)) * fTrans; % nt/sec/ribo
model.vol_Ecoli = 2 ^ ( model.hostGrowth + 2 ) * 1e-16; % liters
model.k_transcriptionEcoli =( 108.7 - 63.6 * model.hostGrowth ^ (-0.1797) );%* fTranscr; % nt/sec/RNAP
model.n_polymeraseEcoli = (192.2 * (model.hostGrowth ) ^ 3 - 155.5 * ( model.hostGrowth ) ^ 2 + 549 * ( model.hostGrowth ) - 103.6 );%* fPoly;
model.n_ribosome = 0.8 * (13.1*model.hostGrowth^2 - 6.89*model.hostGrowth + 6.46) * 1e3 * fRibo;

% rate of genome degradation calculation        
model.k_degradationEcoliDNA = model.NhostGenomePerCell * model.lengthHostGenome * model.fractHostGenomeDegrade * 2 ...
                                / model.tTotalHostGenomeDegrade;
% Note, that in Endy & You this rate is in bp/sec not nt/sec, however for
% interaction with FBA this must be in terms of nt/sec (hence factor of 2)
        
% unit conversions
model.k_packagingT7DNA = model.k_packagingT7DNA / 60 ; % from /min to /sec
model.k_assemblyT7procapsid = model.k_assemblyT7procapsid / 60 ; % from /min to /sec
% model.k_translation = model.k_translation / 3; % from AA/sec/ribosome to nt/sec/ribosome
    
% other parameters
model.tDelayInitial = 0; % sec [code of T7TRANSLOC Endy1997(Thesis)] 
model.switchTranslocation = {'ejectDNAstop' 'phi1o1A'}; % [text Endy1997(Biotech&Bioeng)] // EWB - assume whole element must enter, thus .ntEnd (below) 




% Genome of T7 phage
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% Gene type labels:
% gene = gene
% leftEcRNAPpro = left promotor for E. coli RNA polymerase
% rightEcRNAPpro = right promotor for E. coli RNA polymerase
% termEcRNAP = terminator for E. coli polymerase
% termT7RNAP = terminator for T7 polymerase
% rightT7RNAPpro = right promotor for T7 RNA polymerase
% RNaseIIIsite = site of RNAse activity
% other = other, see notes
%
% [code of T7V2GENOME Endy1997(Thesis)]
% stoichiometry of particle []
% comments from source unless noted with initials EWB
model.lengthT7 = 39937; %nt
model.genomeT7 = struct('A0',     struct( 'type',     'other',... % changed to other since pro strength 0
                                    'ntStart',  224,...
                                    'ntEnd',    289,...
                                    'promotorStrength',	0 ),...
                  'phiOL',  struct( 'type',     'other',... % no transcription
                                    'ntStart',	405,...
                                    'ntEnd',	440 ),...
                  'A1',     struct( 'type',     'rightEcRNAPpro',...
                                    'ntStart',  498,...
                                    'ntEnd',    563,...
                                    'promotorStrength', 1),...
                  'A2',     struct( 'type',     'rightEcRNAPpro',...
                                    'ntStart',  626,...
                                    'ntEnd',    691,...
                                    'promotorStrength', 1),...
                  'A3',     struct( 'type',     'rightEcRNAPpro',...
                                    'ntStart',  750,...
                                    'ntEnd',    815,...
                                    'promotorStrength', 1 ),...
                  'ejectDNAstop',   struct( 'type',     'other',...
                                    'ntStart',  850,...
                                    'ntEnd',    850),...
                  'R0o3',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  816,...
                                    'ntEnd',    890,...
                                    'ntRNase',  890 ),...
              ... % Class I genes
                  'gp0o3',  struct( 'type',     'gene',...
                                    'ntStart',  925,...
                                    'ntEnd',    1278 ),... % 1278, increased to include gp0o4
                  'gp0o4',  struct( 'type',     'gene',...
                                    'ntStart',  1278,...
                                    'ntEnd',    1433 ),...
                  'R0o5',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  1434,...
                                    'ntEnd',    1468,...
                                    'ntRNase',  1468 ),...
                  'gp0o5',  struct( 'type',     'gene',...
                                    'ntStart',  1496,...
                                    'ntEnd',    1639 ),...
                  'B',      struct( 'type',     'other',... % not considered active EcRNAP promotor, in coding region of gp0o5
                                    'ntStart',  1514,...
                                    'ntEnd',    1579 ),...
                  'gp0o6B',  struct( 'type',     'gene',... % part of gp0o5 element
                                    'ntStart',  1636,...
                                    'ntEnd',    1972),...
                  'gp0o6A',  struct( 'type',     'gene',... % part of gp0o5 element
                                    'ntStart',  1636,...
                                    'ntEnd',    1797),...
                  'gp0o7',  struct( 'type',     'gene',... % protien kinase
                                    'ntStart',  2021,...
                                    'ntEnd',    3100),...
                  'C',      struct( 'type',     'other',... % not considered active EcRNAP promotor, in coding region of gp0o5
                                    'ntStart',  3113,...
                                    'ntEnd',    3130 ),...
                  'R1',     struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  3131,... 
                                    'ntEnd',    3138,...
                                    'ntRNase',  3138 ),...
                  'gp1',    struct( 'type',     'gene',... % T7 RNAP
                                    'ntStart',  3171,... 3171,... 
                                    'ntEnd',    5822),...
                  'phi1o1A',struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  5848,...
                                    'ntEnd',    5883,...
                                    'promotorStrength', 0.04),...
                  'R1o1',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  5884,...
                                    'ntEnd',    5887,...
                                    'ntRNase',  5887 ),...
                  'phi1o1B',struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  5923,...
                                    'ntEnd',    5958,...
                                    'promotorStrength', 0.123),...
                  'gp1o1',	struct( 'type',     'gene',... % 
                                    'ntStart',  6007,...
                                    'ntEnd',    6135),...
                  'gp1o2',	struct( 'type',     'gene',... %  inhibits host dGTP triphophohydrolase/target for F plasmid eclusion
                                    'ntStart',  6137,...
                                    'ntEnd',    6394),...
                  'phi1o3',struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  6409,...
                                    'ntEnd',    6444,...
                                    'promotorStrength', 0.007),...
                  'R1o3',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  6445,...
                                    'ntEnd',    6448,...
                                    'ntRNase',  6448 ),...
                  'gp1o3',	struct( 'type',     'gene',... %  DNA ligase
                                    'ntStart',  6475,...
                                    'ntEnd',    7554),...
                  'TE',     struct( 'type',     'termEcRNAP',... 
                                    'ntStart',  7555,...
                                    'ntEnd',    7588,...
                                    'ntTerm',   7588,...
                                    'terminatorStrength', 0.99),...
              ... % Class II genes
                  'gp1o4',	struct( 'type',     'gene',... 
                                    'ntStart',  7608,...
                                    'ntEnd',    7763),...
                  'phi1o5',struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  7778,...
                                    'ntEnd',    7790,...
                                    'promotorStrength', 0.01),...
                  'gp1o5',	struct( 'type',     'gene',... 
                                    'ntStart',  7791,...
                                    'ntEnd',    7880),...
                  'phi1o6',struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  7895,...
                                    'ntEnd',    7905,...
                                    'promotorStrength', 0.01),... 
                  'gp1o6',	struct( 'type',     'gene',... 
                                    'ntStart',  7906,...
                                    'ntEnd',    8166),... % 8166, increased to include gp1o8
                  'gp1o7',	struct( 'type',     'gene',... 
                                    'ntStart',  8166,...
                                    'ntEnd',    8756),...
                  'gp1o8',	struct( 'type',     'gene',... 
                                    'ntStart',  8749,...
                                    'ntEnd',    8895),... % 8895 to eliminate gap
                  'gp2',	struct( 'type',     'gene',... % inhibits E. coli RNAP
                                    'ntStart',  8898,...
                                    'ntEnd',    9092),...
                  'phi2o5', struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  9107,...
                                    'ntEnd',    9142,...
                                    'promotorStrength', 0.01),...
                  'gp2o5',	struct( 'type',     'gene',... % SS DNA binding protein
                                    'ntStart',  9158,...
                                    'ntEnd',    9856),...
                  'gp2o8',	struct( 'type',     'gene',... 
                                    'ntStart',  9857,...
                                    'ntEnd',    10276),... % 10276, extended to include next two
                  'gp3',	struct( 'type',     'gene',...
                                    'ntStart',  10257,...
                                    'ntEnd',    10706),...
                  'gp3o5',	struct( 'type',     'gene',...
                                    'ntStart',  10706,...
                                    'ntEnd',    11161),...
                  'phi3o8', struct( 'type',     'rightT7RNAPpro',...
                                    'ntStart',  11180,...
                                    'ntEnd',    11201,...
                                    'promotorStrength', 0.01),...
                  'R3o8',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  11202,...
                                    'ntEnd',    11203,...
                                    'ntRNase',  11203 ),...
                  'gp3o8',	struct( 'type',     'gene',... % SS DNA binding protein
                                    'ntStart',  11225,...
                                    'ntEnd',    11590),... % 11590, increased to include next 5
                  'gp4a',	struct( 'type',     'gene',... % primase/helicase, part of gp3o8 element
                                    'ntStart',  11565,...
                                    'ntEnd',    13265),...
                  'gp4o1',	struct( 'type',     'gene',... % part of gp3o8 element
                                    'ntStart',  11635,...
                                    'ntEnd',    11757),...
                  'gp4b',	struct( 'type',     'gene',... % helicase, part of gp3o8 element
                                    'ntStart',  11754,...
                                    'ntEnd',    13265),...
                  'phi4c', struct( 'type',     'rightT7RNAPpro',... % part of gp3o8 element
                                    'ntStart',  12671,...
                                    'ntEnd',    12706,...
                                    'promotorStrength', 0.01),...
                  'gp4o2',	struct( 'type',     'gene',... % part of gp3o8 element
                                    'ntStart',  12988,...
                                    'ntEnd',    13326),...   
                  'phi4o3', struct( 'type',     'rightT7RNAPpro',... % part of gp3o8 element
                                    'ntStart',  13341,...
                                    'ntEnd',    13351,...
                                    'promotorStrength', 0.01),...
                  'gp4o3',	struct( 'type',     'gene',... 
                                    'ntStart',  13352,...
                                    'ntEnd',    13564),...
                  'gp4o5',	struct( 'type',     'gene',... % part of gp3o8 element
                                    'ntStart',  13584,...
                                    'ntEnd',    13853),...
                  'R4o7',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  13854,...
                                    'ntEnd',    13892,...
                                    'ntRNase',  13982 ),...
                  'phi4o7', struct( 'type',     'rightT7RNAPpro',... 
                                    'ntStart',  13915,...
                                    'ntEnd',    13926,...
                                    'promotorStrength', 0.01),...
                  'gp4o7',	struct( 'type',     'gene',... 
                                    'ntStart',  13927,...
                                    'ntEnd',    14334),...
                  'gp5',	struct( 'type',     'gene',... 
                                    'ntStart',  14353,...
                                    'ntEnd',    16467),...
                  'gp5o3',	struct( 'type',     'gene',... 
                                    'ntStart',  16483,...
                                    'ntEnd',    16839),...
                  'gp5o5',	struct( 'type',     'gene',... % permits growth on lambda lysogens
                                    'ntStart',  16851,...
                                    'ntEnd',    17359),... % 17150, lengthened to include next 4
                  'gp5o5A',	struct( 'type',     'gene',... % permits growth on lambda lysogens
                                    'ntStart',  16851,...
                                    'ntEnd',    17150),... % 17150, lengthened to include next 4
                  'gp5o7',	struct( 'type',     'gene',... % part of gp5.5 element
                                    'ntStart',  17150,...
                                    'ntEnd',    17359),...
                  'gp5o9',	struct( 'type',     'gene',... % inhibits hose recBCD nuclease, part of gp5.5 element
                                    'ntStart',  17359,...
                                    'ntEnd',    17517),...
                  'gp6',	struct( 'type',     'gene',... % exonuclease, part of gp5.5 element
                                    'ntStart',  17504,...
                                    'ntEnd',    18406),...
                  'gp6o3',	struct( 'type',     'gene',... % part of gp5.5 element
                                    'ntStart',  18394,...
                                    'ntEnd',    18507),...
                  'phi6o5', struct( 'type',     'rightT7RNAPpro',... 
                                    'ntStart',  18545,...
                                    'ntEnd',    18561,...
                                    'promotorStrength', 0.05),...
                  'R6o5',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  18562,...
                                    'ntEnd',    18563,...
                                    'ntRNase',  18563 ),... 
                  'gp6o5',	struct( 'type',     'gene',...  
                                    'ntStart',  18605,...
                                    'ntEnd',    18859),...
                  'gp6o7',	struct( 'type',     'gene',...  
                                    'ntStart',  18864,...
                                    'ntEnd',    19130,...
                                    'stoichPhage', 18),... % 19130, lengthend to include next
                  'gp7',	struct( 'type',     'gene',... % host range
                                    'ntStart',  19130,...
                                    'ntEnd',    19531),... 
                  'gp7o3',	struct( 'type',     'gene',... % host range
                                    'ntStart',  19535,...
                                    'ntEnd',    19834,...
                                    'stoichPhage', 33),...
                  'gp7o7',	struct( 'type',     'gene',... 
                                    'ntStart',  19848,...
                                    'ntEnd',    20240),... % 20240 extended to include next
                  'gp8',	struct( 'type',     'gene',... % head tail connector protien
                                    'ntStart',  20240,...
                                    'ntEnd',    21850,...
                                    'stoichPhage', 12),... 
                  'phi9', struct( 'type',     'rightT7RNAPpro',... 
                                    'ntStart',  21865,...
                                    'ntEnd',    21949,...
                                    'promotorStrength', 0.2),...
                  'gp9',	struct( 'type',     'gene',... % scaffolding protien
                                    'ntStart',  21950,...
                                    'ntEnd',    22873,...
                                    'stoichProcapsid',  137),...
                  'phi10', struct( 'type',     'rightT7RNAPpro',... 
                                    'ntStart',  22904,...
                                    'ntEnd',    22966,...
                                    'promotorStrength', 0.681),...
                  ...'gp10B',	struct( 'type',     'gene',... % minor capsid protien
                  ...                  'ntStart',  22967,...
                  ...                  'ntEnd',    24162),...
                  'gp10A',	struct( 'type',     'gene',... % major capsid protien
                                    'ntStart',  22967,...
                                    'ntEnd',    24004,...
                                    'stoichPhage', 415),... 
                  'TT7',  	struct( 'type',     'termT7RNAP',... 
                                    'ntStart',  24163,...
                                    'ntEnd',    24210,...
                                    'ntTerm',   24210,...
                                    'terminatorStrength', 0.85),...
                  'gp11',	struct( 'type',     'gene',... % tail protien
                                    'ntStart',  24228,...
                                    'ntEnd',    24818,...
                                    'stoichPhage', 12),...
                  'gp12',	struct( 'type',     'gene',... % tail protien
                                    'ntStart',  24842,...
                                    'ntEnd',    27226,...
                                    'stoichPhage', 6),... 
                  'phi13', struct( 'type',     'rightT7RNAPpro',... 
                                    'ntStart',  27274,...
                                    'ntEnd',    27279,...
                                    'promotorStrength', 0.1),...
                  'R13',   struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  27280,...
                                    'ntEnd',    27281,...
                                    'ntRNase',  27281 ),...
                  'gp13',	struct( 'type',     'gene',... % internal virion protien
                                    'ntStart',  27307,...
                                    'ntEnd',    27723),...
                  'gp14',	struct( 'type',     'gene',... % internal virion protien
                                    'ntStart',  27728,...
                                    'ntEnd',    28318,...
                                    'stoichPhage', 10),...
                  'gp15',	struct( 'type',     'gene',... % internal virion protien
                                    'ntStart',  28325,...
                                    'ntEnd',    30568,...
                                    'stoichPhage', 8),...
                  'gp16',	struct( 'type',     'gene',... % internal virion protien
                                    'ntStart',  30595,...
                                    'ntEnd',    34551,...
                                    'stoichPhage', 4),...
                  'phi17', struct( 'type',     'rightT7RNAPpro',... 
                                    'ntStart',  34566,...
                                    'ntEnd',    34623,...
                                    'promotorStrength', 0.1),...
                  'gp17',	struct( 'type',     'gene',... % tail fiber protein
                                    'ntStart',  34624,...
                                    'ntEnd',    36285,...
                                    'stoichPhage', 18),...
                  'gp17o5',	struct( 'type',     'gene',... % lysis protein
                                    'ntStart',  36344,...
                                    'ntEnd',    36547),...
                  'gp18',	struct( 'type',     'gene',... % DNA maturation protein
                                    'ntStart',  36553,...
                                    'ntEnd',    36822),...
                  'R18o5', 	struct( 'type',     'RNaseIIIsite',...
                                    'ntStart',  36855,...
                                    'ntEnd',    36856,...
                                    'ntRNase',  36856 ),...
                  'gp18o5',	struct( 'type',     'gene',... % homologous to lambda Rz lysis protein
                                    'ntStart',  36917,...
                                    'ntEnd',    37348),...
                  'gp18o7',	struct( 'type',     'gene',... %  within coding region of 18o5
                                    'ntStart',  37032,...
                                    'ntEnd',    37283),...
                  'gp19',	struct( 'type',     'gene',... % DNA maturation protein
                                    'ntStart',  37370,...
                                    'ntEnd',    39130),...
                  'gp19o2',	struct( 'type',     'gene',... % within coding region of 19
                                    'ntStart',  38016,...
                                    'ntEnd',    38273),...
                  'gp19o3',	struct( 'type',     'gene',... % within coding region of 19
                                    'ntStart',  38553,...
                                    'ntEnd',    38726),...
                  'phiOR', struct(  'type',     'other',... 
                                    'ntStart',  39229,...
                                    'ntEnd',    39388),...
                  'gp19o5',	struct( 'type',     'gene',... % within coding region of 19
                                    'ntStart',  39389,...
                                    'ntEnd',    39538)); % nt
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 



% Set maximum polymerase density
% model.maxDensityRNApolymerase = 1 / 200; % as per Thesis code

% Genome translocation, for calculation of X
% start time of translocation mediated by E. coli RNA polymerase
model.tStartEcoliRNAP = model.tDelayInitial +  model.genomeT7.(model.switchTranslocation{1}).ntEnd ...
                        / model.k_translocationT7(1) ; % sec                 
% start time of translocation mediated by T7 RNA polymerase                   
model.tStartT7RNAP = model.tStartEcoliRNAP + ( model.genomeT7.(model.switchTranslocation{2}).ntEnd - ...
                 model.genomeT7.(model.switchTranslocation{1}).ntEnd ) / model.k_translocationT7(2); % sec
% translocation completion time
model.tFinishTranslocation = model.tStartT7RNAP +  ( model.lengthT7 - ...
                 model.genomeT7.(model.switchTranslocation{2}).ntEnd ) / model.k_translocationT7(3); % sec


            
%% Collect genetic elements by type
% This sorts the genetic elements by their type as specified in the
% section above "Genome T7 Phage".
model = sortGeneticElementsByType(model);
%% Collect genes depending on which mRNA they will appear on
% %% Collect genes depending on which mRNA they will appear on
[mRNA,gene, mRNAcompact] = mRNAgeneratorIndividualmRNA(model);
% [mRNA,gene] = mRNAgenerator(model);
model.mRNA = mRNA;
model.n_mRNA = size(mRNA,1);
model.gene = gene;
model.mRNAcompact = mRNAcompact;

% %% Plots possible mRNA's with genes on them
% plotsT7Genome(model);

%% some calculations for setup

% 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
model.downstreamOfTermT7RNAPpro = (model.rightT7RNAPpro{3} > model.genomeT7.TT7.ntTerm );


model.mRNAupstreamT7RNAPterm = false( model.n_gene,1);
model.mRNAdownstreamT7RNAPterm = false( model.n_gene,1);
model.mRNAupstreamEcRNAPterm = false( model.n_gene,1);
for k = 1 : model.n_gene
    
    % 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)
    model.mRNAupstreamT7RNAPterm(k) =  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)
    model.mRNAdownstreamT7RNAPterm(k) =  model.genomeT7.TT7.ntTerm < 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)
    model.mRNAupstreamEcRNAPterm(k) =  model.genomeT7.TE.ntTerm >= model.mRNA{k,4};

    % 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
    model.mRNAdownstreamEcRNAPpro{k} = (model.rightEcRNAPpro{2}(:,1)) <= model.mRNA{k,3};
    
    % 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
    model.mRNAdownstreamT7RNAPpro{k} = (model.rightT7RNAPpro{2}(:,1)) <= model.mRNA{k,3};
    
    
end



end % main function


% subfunctions
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% sortGeneticElementsByType - Sorts the elements of T7 genome based on
% their "type" as specified below
%
% Programmer:
%   Elsa Birch - ewbirch@stanford.edu
%
% Gene type labels:
% gene = gene
% leftEcRNAPpro = left promotor for E. coli RNA polymerase
% rightEcRNAPpro = right promotor for E. coli RNA polymerase
% termEcRNAP = terminator for E. coli polymerase
% termT7RNAP = terminator for T7 polymerase
% rightT7RNAPpro = right promotor for T7 RNA polymerase
% RNaseIIIsite = site of RNAse activity
% other = other, see notes

function [model] = sortGeneticElementsByType(model)

geneticElements = fieldnames(model.genomeT7);
n_geneticElements = length(geneticElements);

% initialize storage cell arrays:
% NOTE: sizes not known beforehand, resize durring sort, could be
% preallocated but only occurs once
gene = {}; % name (str), start (nt), end (nt), genome stoich (#), particle stoich (#)
leftEcRNAPpro = {};
rightEcRNAPpro = {};
termEcRNAP = {};
termT7RNAP = {};
rightT7RNAPpro = {};
RNaseIIIsite = {};
geneAminoAcids = [];

% counters for each type of element
n_gene = 0; 
n_leftEcRNAPpro = 0;
n_rightEcRNAPpro = 0;
n_termEcRNAP = 0;
n_termT7RNAP = 0;
n_rightT7RNAPpro = 0;
n_RNaseIIIsite = 0; 

% loop through genetic elements and sort based on type field
for k = 1:n_geneticElements
    switch model.genomeT7.(geneticElements{k}).type
        case 'gene'
            n_gene = n_gene +1; % update counter
            gene{1}(n_gene,1) = {geneticElements{k}};
            gene{2}(n_gene,1) = model.genomeT7.(geneticElements{k}).ntStart;
            gene{3}(n_gene,1) = model.genomeT7.(geneticElements{k}).ntEnd;
            gene{4}(n_gene,1) = 0; % assume not components of particle, correct below
            gene{5}(n_gene,1) = 0;
            % if GP is part of PC or assembled phage, store stoich
            % coefficient
            if isfield(model.genomeT7.(geneticElements{k}),'stoichPhage')
                gene{4}(n_gene,1) = model.genomeT7.(geneticElements{k}).stoichPhage;
            elseif isfield(model.genomeT7.(geneticElements{k}),'stoichProcapsid')
                gene{5}(n_gene,1) = model.genomeT7.(geneticElements{k}).stoichProcapsid;
            end %if
            
            % store the location in gene cell array in model structure for
            % easy gene name --> property reference
            model.genomeT7.(geneticElements{k}).n_gene = n_gene;

        case 'leftEcRNAPpro'
            n_leftEcRNAPpro = n_leftEcRNAPpro + 1;
            leftEcRNAPpro(n_leftEcRNAPpro,4) = {geneticElements{k},...
                                model.genomeT7.(geneticElements{k}).ntStart,...
                                model.genomeT7.(geneticElements{k}).ntEnd,...
                                model.genomeT7.(geneticElements{k}).promotorStrength};
        case 'rightEcRNAPpro'
            n_rightEcRNAPpro = n_rightEcRNAPpro + 1;
            rightEcRNAPpro{1}{n_rightEcRNAPpro,1} = geneticElements{k};
            rightEcRNAPpro{2}(n_rightEcRNAPpro,1) = model.genomeT7.(geneticElements{k}).ntStart;
            rightEcRNAPpro{3}(n_rightEcRNAPpro,1) = model.genomeT7.(geneticElements{k}).ntEnd;
            rightEcRNAPpro{4}(n_rightEcRNAPpro,1) = model.genomeT7.(geneticElements{k}).promotorStrength ;
%             rightEcRNAPpro(n_rightEcRNAPpro,1:4) = {geneticElements{k},...
%                                 model.genomeT7.(geneticElements{k}).ntStart,...
%                                 model.genomeT7.(geneticElements{k}).ntEnd,...
%                                 model.genomeT7.(geneticElements{k}).promotorStrength};
        case 'termEcRNAP'
            n_termEcRNAP = n_termEcRNAP + 1;
            termEcRNAP(n_termEcRNAP,1:5) = {geneticElements{k},...
                                model.genomeT7.(geneticElements{k}).ntStart,...
                                model.genomeT7.(geneticElements{k}).ntEnd,...
                                model.genomeT7.(geneticElements{k}).ntTerm,...
                                model.genomeT7.(geneticElements{k}).terminatorStrength};
        case 'termT7RNAP'
            n_termT7RNAP = n_termT7RNAP + 1;
            termT7RNAP(n_termT7RNAP,1:5) = {geneticElements{k},...
                                model.genomeT7.(geneticElements{k}).ntStart,...
                                model.genomeT7.(geneticElements{k}).ntEnd,...
                                model.genomeT7.(geneticElements{k}).ntTerm,...
                                model.genomeT7.(geneticElements{k}).terminatorStrength};
        case 'rightT7RNAPpro'
            n_rightT7RNAPpro = n_rightT7RNAPpro + 1;
            rightT7RNAPpro{1}{n_rightT7RNAPpro,1} = geneticElements{k};
            rightT7RNAPpro{2}(n_rightT7RNAPpro,1) = model.genomeT7.(geneticElements{k}).ntStart;
            rightT7RNAPpro{3}(n_rightT7RNAPpro,1) = model.genomeT7.(geneticElements{k}).ntEnd;
            rightT7RNAPpro{4}(n_rightT7RNAPpro,1) = model.genomeT7.(geneticElements{k}).promotorStrength ;
%             rightT7RNAPpro(n_rightT7RNAPpro,1:4) = {geneticElements{k},...
%                                 model.genomeT7.(geneticElements{k}).ntStart,...
%                                 model.genomeT7.(geneticElements{k}).ntEnd,...
%                                 model.genomeT7.(geneticElements{k}).promotorStrength};
        case 'RNaseIIIsite'
            n_RNaseIIIsite = n_RNaseIIIsite + 1;
            RNaseIIIsite(n_RNaseIIIsite,1:4) = {geneticElements{k},...
                                model.genomeT7.(geneticElements{k}).ntStart,...
                                model.genomeT7.(geneticElements{k}).ntEnd,...
                                model.genomeT7.(geneticElements{k}).ntRNase};
        case 'other'
            % do nothing
        otherwise
            % warn if type is unrecognized
            disp(['unknown genetic element type: ', ...
                 model.genomeT7.(geneticElements{k}).type])
    end % switch
end % for k

% Store results/counts in model structure
model.gene = gene;
model.leftEcRNAPpro = leftEcRNAPpro;
model.rightEcRNAPpro = rightEcRNAPpro;
model.termEcRNAP = termEcRNAP;
model.termT7RNAP = termT7RNAP;
model.rightT7RNAPpro = rightT7RNAPpro;
model.RNaseIIIsite = RNaseIIIsite;
model.n_gene = n_gene;
model.n_leftEcRNAPpro = n_leftEcRNAPpro;
model.n_rightEcRNAPpro = n_rightEcRNAPpro;
model.n_termEcRNAP = n_termEcRNAP;
model.n_termT7RNAP = n_termT7RNAP;
model.n_rightT7RNAPpro = n_rightT7RNAPpro;
model.n_RNaseIIIsite = n_RNaseIIIsite;
model.geneAminoAcids = geneAminoAcids;


end % function sortGeneticElementsByType


function [mRNAfinal,gene,mRNAcompact] = mRNAgeneratorIndividualmRNA(model)

% preallocate memory for all possible mRNA, more than needed, but truncate
% later to avoid changeing variable size on the fly
n_mRNAallCombo = 2 * (model.n_RNaseIIIsite + model.n_rightT7RNAPpro + ...
    model.n_termEcRNAP + model.n_termT7RNAP + model.n_rightEcRNAPpro);
mRNA{n_mRNAallCombo-1,4} = [];

mRNAcount = 0;
geneticElements = fieldnames(model.genomeT7);

for k = 1:length(geneticElements)
    if strcmp(model.genomeT7.(geneticElements{k}).type,'gene')
        mRNAcount = mRNAcount + 1;
        mRNA(mRNAcount,1:2) = {geneticElements{k},geneticElements{k}};
        mRNA(mRNAcount,3:4) = {model.genomeT7.(geneticElements{k}).ntStart, model.genomeT7.(geneticElements{k}).ntEnd};
    end
end
mRNAfinal = mRNA;
% Add rest of information on mRNA
model.gene{6} = [];
for iGene = 1:model.n_gene
    model.gene{6}(iGene,1) = iGene;
end

gene = model.gene;

% create 1x4 cell array of column vectors/cells
mRNAcompact = { {mRNA{:,1}}', {mRNA{:,2}}',  cell2mat({mRNA{:,3}})', cell2mat({mRNA{:,4}})'};

end % function mRNAgeneratorIndividualmRNA