% Author: Elsa W. Birch, ewbirch@stanford.edu, Covert Lab - Stanford Bioengineering
% Date: 4/24/2011
% Version: associated with production of the figures for Birch et al. 2012 
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Summary: parses organism to simulation structure
% 
%
% Inputs: below, (may be out of date)chck if up to date)
%
% Outputs: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function [ organismModelStruct, organismModelWarnings] = ...
                    FBA_EB_organismFileParse( organismModelFolder, organismModelName, organismParseOptions )
% inputs - organismModelFolder: Full path to file in folder containing organism model,
%                               including the name of the folder itself.
%                               Naming of folder should be 'nameModel', and should contain 
%                               tab delimited files: nameReactions.txt
%                                                    nameMetabolites.txt
%                                                    nameBiomass.txt
%        - organismParseOptions: Structure of flags/options for parsing:
%                                       .fluxBoundsDefault is a two element array
%                                        containing the [lower upper] bounds to be used if
%                                        none other are set for a given reaction
%                                       .geneAutoAddOptions.maxGeneNumber
%                                       .geneAutoAddOptions.flag 
%                                       .geneAutoAddOptions.delayDefault
%                                       .regulationFlag, true ( 1 ) if regualtion used,
%                                        false ( 0 ) if regulation will be ignored
%
% ( format of files expected in organismModelFolder ) 
%   - nameReactions.txt: tab deleimited text, with columns [reaction name], [reaction
%                        equation], [comments]. The name of reactions should include
%                        alphanumerice characters and underscores only, and should be
%                        unique (dashes will be replaced with underscores, commas will be
%                        removed).  Reactions should be preceeded by a compartment label 
%                        contained in square brackets, followed by a space colon space. 
%                        Reactions should include reactants preceeded by a
%                        stoichiometric coefficient and space, each separated by a + sign,
%                        then products after a separator, either ==> for
%                        irreversible
%                        reactions or <==> for reversible reactions. Example:
%                         [c] : 1 A + 2 B ==> 1 D  
%                        alternate example:
%                        
%
%
%
% Conventions for naming etc
%       Compartments, the designation that is automatically assigned special reactions
%                     should be avoided unless intended:
%                     t - transport reactions (those more than 1 compartment). This is not
%                         a designation that is applied to metabolites.
%                     e - environment, extracellular, export, metabolites. Reactions
%                         associated with these metabolites are typically transport.
%                     a - 
%

% Note that the structure of error collection and output is designed to output the most
% possible warning/error information for each parsing run. This is intended to be useful,
% such that (for example) multiple equation entry errors are identified and output in one
% run, rather than requireing a run to identify every error.
%

% Initial preperation for and load of files into raw text format
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% file prefix is letters preceeding 
organismModelNamePrefix = organismModelName; 
% Reactions File
fileIdRxns = fopen( fullfile( organismModelFolder, [organismModelNamePrefix 'Reactions.txt']) );
% read in contents, assuming format columns: rxnName, rxnDefinition, rxnRegulation,
% metabolic classification, metabolic sub-classification, 4x(comment(discard))
headerRowsRxns = 1; % # ignored header rows
rawTextReadRxns = textscan(fileIdRxns, '%s %s %s %s %s %s %*s %*s %*s %*s', 'Delimiter', '\t','HeaderLines', headerRowsRxns );
fclose(fileIdRxns); % close file

% Biomass file
fileIdBiomass = fopen(  fullfile( organismModelFolder, [organismModelNamePrefix 'Biomass.txt'] ));
% read in contents, assuming format columns: metaboliteName, fluxUpper, fluxLower, comment(discard)
headerRowsBiomass = 1; % # ignored header rows
rawTextReadBiomass = textscan(fileIdBiomass, '%s %s %f %*s', 'Delimiter', '\t', 'HeaderLines', headerRowsBiomass );
fclose(fileIdBiomass); % close file

% Flux bounds file
fileIdFluxBounds = fopen( fullfile( organismModelFolder, [organismModelNamePrefix 'FluxBounds.txt'] ));
% read in contents, assuming format columns: rxnName, fluxUpper, fluxLower, fulxCondition, comment(discard)
headerRowsFluxBounds = 1; % # ignored header rows
rawTextReadFluxBounds = textscan(fileIdFluxBounds, '%s %s %f %f %s %*s %*s %*s %*s', 'Delimiter', '\t','HeaderLines', headerRowsFluxBounds );
fclose(fileIdFluxBounds); % close file

% Regulatory files
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

% Gene file
fileIdGenes = fopen( fullfile( organismModelFolder,  [organismModelNamePrefix 'Genes.txt']) );
% read in contents, assuming format columns: geneNumber, geneName, regulationRule, comment(discard)
headerRowsGenes = 1; % # ignored header rows
rawTextReadGenes = textscan(fileIdGenes, '%s %s %s %f %*s %*s %*s %*s', 'Delimiter', '\t','HeaderLines', headerRowsGenes );
fclose(fileIdGenes); % close file

% load from TF file
fileIdTFs = fopen( fullfile( organismModelFolder, [organismModelNamePrefix 'TFs.txt']) );
% read in contents, assuming format columns: nameTF, TFregulationRule, comment(discard)
headerRowsTFs = 1; % # ignored header rows
rawTextReadTFs = textscan(fileIdTFs, '%s %s %*s %*s %*s %*s', 'Delimiter', '\t', 'HeaderLines', headerRowsTFs );
fclose(fileIdTFs); % close file

% load from Stimuli file
fileIdStimuli = fopen( fullfile( organismModelFolder,  [organismModelNamePrefix 'Stimuli.txt']) );
% read in contents, assuming format columns: nameStimuli, Stimuli values, comment(discard)
headerRowsStimuli = 1; % # ignored header rows
rawTextReadStimuli = textscan(fileIdStimuli, '%s %s %*s', 'Delimiter', '\t', 'HeaderLines', headerRowsStimuli );
fclose(fileIdStimuli); % close file

% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

% load from Metabolite file
fileIdMetabolites = fopen(fullfile( organismModelFolder,  [organismModelNamePrefix 'Metabolites.txt']) );
% read in contents, assuming format columns: metaboliteAbreviation(model), metaboliteFullName, comment(discard)
headerRowsMetabolites = 1; % # ignored header rows
rawTextReadMetabolites = textscan(fileIdMetabolites, '%s %s %*s', 'Delimiter', '\t', 'HeaderLines', headerRowsMetabolites );
fclose(fileIdMetabolites); % close file


% load from Cycles file
fileIdCycles = fopen(fullfile( organismModelFolder,  [organismModelNamePrefix 'Cycles.txt']) );
% read in contents, assuming format columns: metaboliteAbreviation(model), metaboliteFullName, comment(discard)
headerRowsCycles = 1; % # ignored header rows
rawTextReadCycles= textscan(fileIdCycles, '%s %*s', 'Delimiter', '\t', 'HeaderLines', headerRowsCycles );
fclose(fileIdCycles); % close file

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% ( Initial preperation for and load of files into raw text format )

% initialize structures
organismModelStruct.metabolites.dummy = []; % not ideal way to do this, fix later (.dummy field removed from structure later)
organismModelStruct.metaboliteList = {};
organismModelStruct.nMetabolites = [];
organismModelStruct.rxns.dummy = []; % not ideal way to do this, fix later (.dummy field removed from structure later)
organismModelStruct.rxnList = {};
organismModelStruct.nRxns = [];
organismModelStruct.stoichMatrix = [];
organismModelStruct.biomass = [];
organismModelStruct.fluxBoundsDefault = [];
organismModelStruct.metaboliteCompartmentList = {};
organismModelStruct.rxnCompartmentList = {};
organismModelStruct.reversibility = [];
organismModelWarnings = {};
organismModelStruct.nStimuli = [];
organismModelStruct.nMediaExchange = [];



% Parse/Construct Basic Reaction & Metabolite Functional Components of the FBA file/model
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% Parse Reactions
% using subfunction (see end of file)
[organismModelStruct, organismModelErrors] = parseRxns(rawTextReadRxns, organismModelStruct);
% throw error and exit if errors have collected so far, because later errors may be
% resultant from these initial problems and therefore comfusing
if ~isempty(organismModelErrors)
    for iEr = 1 : length(organismModelErrors)
            disp(sprintf('%s', organismModelErrors{iEr} ))
    end % iEr
    error( ['failure in parsing', organismModelNamePrefix ,...
            'Reactions.txt file, see above error massages for more explaination' ])
end % throw error if any have been found so far


% Parse Biomass
% uses subfunction, (see end of this file)
[organismModelStruct,organismModelErrors] = parseBiomass(rawTextReadBiomass,...
                                                         organismModelStruct,...
                                                         organismModelErrors);
% throw error and exit if errors have collected so far, because later errors may be
% resultant from these initial problems and therefore comfusing
if ~isempty(organismModelErrors)
    disp( organismModelErrors )
    error( ['failure in parsing ', organismModelNamePrefix ,...
            'Biomass.txt file, see above error massages for more explaination' ])
end % throw error if any have been found so far


% Add Exchange reactions to exchange 'e' metabolites
organismModelStruct = addExchangeReactions(organismModelStruct);


% Add Accumulation reactions to accumulate 'a' metabolites
organismModelStruct = addAccumulationReactions(organismModelStruct);

% Find cycles
[ organismModelStruct ] = findSimpleCycles(organismModelStruct);

% Check for 'orphan' metabolites, and store met -> rxn asscociations
[organismModelStruct, orphanWarnings] = checkMetaboliteOrphans(organismModelStruct);
% throw error and exit if errors have collected so far, because later errors may be
% resultant from these initial problems and therefore comfusing
if ~isempty(orphanWarnings)
    disp(['WARNING: in check of stoichiometric matrix feasibility, see below massages for explaination'])
    for iWarn = 1 : length(orphanWarnings)
        disp(sprintf('%s',  orphanWarnings{iWarn} ))
    end % iWarn
    disp('WARNING: in check of stoichiometric matrix feasibility, see above massages for explaination')
    disp(' ')
end % throw error if any have been found so far


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% ( Parse Basic Functional Components of the FBA file )







    
    


% Parse/Construct regulatory portions of model
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if organismParseOptions.regulationFlag

    % Parse Genes 
    [organismModelStruct, geneErrors] = ...
        parseGenes(rawTextReadGenes, organismModelStruct, organismParseOptions.geneAutoAddOptions);

    % Parse TFs
    [organismModelStruct, tfErrors] = parseTFs(rawTextReadTFs, organismModelStruct);

    % Parse Stimuli
    [organismModelStruct, stimuliErrors] = parseStimuli(rawTextReadStimuli, organismModelStruct);

    % throw error and exit if errors have collected so far, because later errors may be
    % resultant from these initial problems and therefore comfusing
    organismModelErrors = [organismModelErrors; geneErrors; tfErrors; stimuliErrors];
    if ~isempty(organismModelErrors)
        organismModelErrors
        error( 'organism setup failure, see above error massages for more explaination' )
    end % throw error if any have been found so far

    % Parse TF rules
    organismModelStruct.TFgeneActivityInstant = organismParseOptions.TFgeneActivityInstant;
    [organismModelStruct, tfRuleErrors] = parseTFsRules(organismModelStruct);
    organismModelStruct.TFgeneActivityInstant = false;
    
    % Parse Gene rules
    [organismModelStruct, geneRuleErrors] = parseGenesRules(organismModelStruct);

    % Parse Rxn rules
    [organismModelStruct, rxnRuleErrors] = parseRxnsRules(organismModelStruct);
    
    % 
    [organismModelStruct, geneDependErrors] = parseRxnsGenes(organismModelStruct);
    
    % Throw any errors
    organismModelErrors = [organismModelErrors; geneRuleErrors; tfRuleErrors; rxnRuleErrors, geneDependErrors];
    if ~isempty(organismModelErrors)
        for iEr = 1 : length(organismModelErrors)
            disp(sprintf('%s', organismModelErrors{iEr} ))
        end % iEr
        error( 'organism regulation setup failure, see above error massages for more explaination' )
    end % throw error if any have been found so far

end % if organismParseOptions.regualtionFlag

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% (Parse/Construct regulatory portions of model)

% rawTextReadFluxBounds{1}
% rawTextReadFluxBounds{3}

% Parse Flux Bounds
[organismModelStruct, organismModelErrors] = ...
    parseFluxBounds(rawTextReadFluxBounds, organismModelStruct, organismParseOptions.fluxBoundsDefault);
% throw error and exit if errors have collected so far, because later errors may be
% resultant from these initial problems and therefore comfusing
if ~isempty(organismModelErrors)
    for iEr = 1 : length(organismModelErrors)
        disp(sprintf('%s', organismModelErrors{iEr} ))
    end % iEr
    error( 'failure in flux bound parsing, see above for detail' )
end % throw error if any have been found so far

% parse cycles file
[organismModelStruct] = parseCycles( organismModelStruct ,rawTextReadCycles );

% Labeling Parsing 
% Note: uneccesary for basic functionality of model, however may be refferenced in other 
% code, thus only warnings (not errors) are produced in this section
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 


[organismModelStruct, unnamedMet] = ...
    parseMetaboliteNames(rawTextReadMetabolites, organismModelStruct);
% disp( sprintf( '%s', 'the following metabolite abbreviations lack an associated full chemical name'))
% disp( unnamedMet )


[organismModelStruct ] = parseMetabolicClass( organismModelStruct );

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

end % function FBA_EB_organismFileParse


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

function [organismModelStruct, organismModelErrors] = parseRxns(rawTextReadRxns, organismModelStruct)

% initialize error array
organismModelErrors = {};

rawTextReadRxnNames = rawTextReadRxns{ 1 };
rawTextReadRxnEqn = rawTextReadRxns{ 2 };
rawTextReadRxnRegulation = rawTextReadRxns{ 3 };
rawTextReadRxnGeneDependance = rawTextReadRxns{ 4 };
rawTextReadRxnMetabolicClass = rawTextReadRxns{ 5 };
rawTextReadRxnMetabolicSubClass = rawTextReadRxns{ 6 };

% catch for uneven column lengths if no rule for last entry
if size(rawTextReadRxnEqn,1) - size(rawTextReadRxnRegulation,1) == 1;
    % one less entry found for reaction rules than reaction name/eqns, assume this was an
    % empty cell not captured in tab delimiting of text and add a whitespace string there
    rawTextReadRxnRegulation = [rawTextReadRxnRegulation; ' '];
end 


% preliminary count of reactions
nRxnPrelim = size( rawTextReadRxnNames , 1);

% initialize valid reaction count, and metabolite count
nRxns = 0;
nMetabolites = 0;

% loop through all reactions loaded, and store
for iRawRxns = 1 : nRxnPrelim
    
    % reset failure flag for adding reaction to model
    rxnFail = 0;
    
    % reaction name, and cleaning
    rawRxnName = rawTextReadRxnNames{ iRawRxns };
    cleanRxnName = strrep(rawRxnName,'-','_'); % replace dash with underscore
    cleanRxnName = strrep(cleanRxnName,',',''); % remove commas
    rxnName = cleanRxnName;
    
    % get and parse reaction
    rawRxnEqn = rawTextReadRxnEqn{ iRawRxns };
    [reactants, products, reversible, compartment, rxnEqnParseFail ] = FBA_EB_rxnEqnParse(rawRxnEqn); 

    % add compartment (as lowercase) to front of rxn name (uppercase), 
    rxnNameWithCompartment = [lower(compartment),upper(rxnName)];
    
    % Tests for valididty and uniqueness of varible names etc, and adds warning as needed
    % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    
    % test rxn name valid for variable name use
    % note: later a lowercase letter designating compartment will preceede this, but still
    % check if reaction name alone is valid for varible use
    if ~isvarname(rxnName)
        % invalid, fail and add warning about this variable name
        rxnFail = 1;
        organismModelErrors = [ organismModelErrors ;
                                 ['organism reaction name: ', rawRxnName,' changed to: ',...
                                 rxnName,' but is still invalid variable name'] ] ;                       
    end % if invalid rxn name
    
    if FBA_EB_reservedKeywordCheck(rxnName)
        % invalid, fail and add warning about this variable name
        rxnFail = 1;
        organismModelErrors = [ organismModelErrors ;
                                 ['organism reaction name: ', rawRxnName,' changed to: ',...
                                 rxnName,' but is still invalid variable name'] ] ;                       
    end % if reserved keyword in rxn name
    
    % test rxn name for uniqueness (against all current)
    if isfield(organismModelStruct.rxns, rxnNameWithCompartment)
        % rxn name not unique, fail and add warning
        rxnFail = 1;
        organismModelErrors = [ organismModelErrors ;
                           ['reaction name is not unique: ', rxnNameWithCompartment ] ] ;
    end % rxn name is not unique
    
    % if failure output from equation parser, add warning
    if rxnEqnParseFail
        rxnFail = 1;
        organismModelErrors = [ organismModelErrors ;
                               ['failure parsing reaction eqn(',rxnNameWithCompartment,'): ', rawRxnEqn ] ] ;
    end % if rxn eqn parse failure
    % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    
    
    % no failure so far, add to stoich matrix
    if ~rxnFail
        
        % add to valid reaction count
        nRxns = nRxns + 1;
        
        % update rxn list and structure in organism structure
        organismModelStruct.rxnList{nRxns} = rxnNameWithCompartment;
        organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn = nRxns; % at some point check that these are all used
        organismModelStruct.rxns.(rxnNameWithCompartment).reactants = reactants; 
        organismModelStruct.rxns.(rxnNameWithCompartment).products = products; 
        organismModelStruct.rxns.(rxnNameWithCompartment).reversible = reversible; 
        organismModelStruct.rxns.(rxnNameWithCompartment).rawRxnString = rawRxnEqn ;
        organismModelStruct.rxns.(rxnNameWithCompartment).rawRegulationString = rawTextReadRxnRegulation{ iRawRxns } ;
        organismModelStruct.rxns.(rxnNameWithCompartment).metabolicClass = rawTextReadRxnMetabolicClass{ iRawRxns };
        organismModelStruct.rxns.(rxnNameWithCompartment).metabolicSubClass = rawTextReadRxnMetabolicSubClass{ iRawRxns };
        organismModelStruct.rxns.(rxnNameWithCompartment).rawGeneDependanceString = rawTextReadRxnGeneDependance{ iRawRxns };
        % store reversibility
        organismModelStruct.reversibility(nRxns) = reversible;
        
        % store compartment
        organismModelStruct.rxnCompartmentList{nRxns} = compartment;

        
        % combine reactants and products to loop through metabolites
        rxnMetabolites = [reactants ; products];
        nRxnMetabolites = size( rxnMetabolites, 1);
        
        % loop through metabolites adding to stoichiometric matrix, and lists as needed
        for iRxnMet = 1 : nRxnMetabolites;
            
            % add compartment to metabolite name
            metaboliteNameWithCompartment = [lower(rxnMetabolites{iRxnMet,2}),...
                                                rxnMetabolites{iRxnMet,1}];
                                            
            % determine if metabolite already exists in model, and add if not
            if ~isfield(organismModelStruct.metabolites, metaboliteNameWithCompartment)
                % metabolite is not yet included in model
                
                % add to metabolite count
                nMetabolites = nMetabolites + 1;
                
                % add to metabolite structure and list
                organismModelStruct.metaboliteList{nMetabolites} = metaboliteNameWithCompartment;
                organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite = nMetabolites ;
                organismModelStruct.metaboliteCompartmentList{nMetabolites} = metaboliteNameWithCompartment(1);
                
            end % if this is a new metabolite
            
            % add metabolite to stoichiometric matrix
            rowStoichMatrix = organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite;
            organismModelStruct.stoichMatrix(rowStoichMatrix, nRxns) = rxnMetabolites{iRxnMet,3};
            
        end % for iRxnMet
        
    end % if no failure add to stoich matrix
    
end % for iRawRxns

% store reaction and metabolite counts
organismModelStruct.nRxns = nRxns;
organismModelStruct.nMetabolites = nMetabolites;

% remove dummy fields from initialization of metabolites and rxns structures
organismModelStruct.metabolites = rmfield(organismModelStruct.metabolites, 'dummy');
organismModelStruct.rxns = rmfield(organismModelStruct.rxns, 'dummy');

end % function parseRxns

function [organismModelStruct,organismModelErrors] = parseBiomass(rawTextReadBiomass,organismModelStruct,organismModelErrors)

organismModelErrors = {};

rawTextReadBiomassMetabolites = rawTextReadBiomass{ 1 };
rawTextReadBiomassCompartments = rawTextReadBiomass{ 2 };
rawTextReadBiomassStoich = rawTextReadBiomass{ 3 };

% initialize biomass metabolite count and storage
nBiomassMetabolites = size(rawTextReadBiomassMetabolites, 1 );
biomassMetabolites{nBiomassMetabolites, 3} = [];

% add one to overall reaction count
organismModelStruct.nRxns = organismModelStruct.nRxns + 1;

% loop through biomass components matabolite adding to reaction
for iBiomassMet = 1 : nBiomassMetabolites
    
    % add compartment (as lowercase) to front of rxn name (uppercase), 
    metaboliteNameWithCompartment = [ lower(rawTextReadBiomassCompartments{ iBiomassMet }), ...
                                        upper(strrep(rawTextReadBiomassMetabolites{ iBiomassMet },'-','_')) ];
    metaboliteStoichCoeff = rawTextReadBiomassStoich(iBiomassMet);
                                    
    % check that this metabolite exists
    if isfield(organismModelStruct.metabolites, metaboliteNameWithCompartment)
        % metabolite exists in model, and can therefore be added as a biomass component
        
        % store for reaction reffernce
        biomassMetabolites{iBiomassMet, 1} = metaboliteNameWithCompartment(2:end);
        biomassMetabolites{iBiomassMet, 2} = metaboliteNameWithCompartment(1);
        biomassMetabolites{iBiomassMet, 3} = metaboliteStoichCoeff;
        
        % add to stoichometric matrix
        rowStoichMatrix = organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite;
        organismModelStruct.stoichMatrix(rowStoichMatrix, organismModelStruct.nRxns) = metaboliteStoichCoeff;
        
    else
        % metabolite not found in model
        organismModelErrors = [ organismModelErrors ;
                               ['biomass metabolite: ', metaboliteNameWithCompartment,' not found in any other rxn' ] ] ;
    
    end % if biomass metabolite is found in model metabolites
    
end % for mBiomassMet

% add biomass metabolite to count, structure, and list
metaboliteNameWithCompartment = 'aBIOMASS';
organismModelStruct.nMetabolites = organismModelStruct.nMetabolites + 1;
organismModelStruct.metaboliteList{organismModelStruct.nMetabolites} = metaboliteNameWithCompartment;
organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite = organismModelStruct.nMetabolites ;
organismModelStruct.metaboliteCompartmentList{organismModelStruct.nMetabolites} = metaboliteNameWithCompartment(1);


% add biomass reaction to stoich matrix
rowStoichMatrix = organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite;
organismModelStruct.stoichMatrix(rowStoichMatrix, organismModelStruct.nRxns) = 1;

% update rxn list and structure in organism structure
rxnNameWithCompartment = 'bBIOMASS';
idxBiomassReactants = cell2mat(biomassMetabolites(:, 3)) < 0 ;
reactants = biomassMetabolites(idxBiomassReactants, :);
idxBiomassProducts = cell2mat(biomassMetabolites(:, 3)) > 0 ;
products = biomassMetabolites(idxBiomassProducts, :);
organismModelStruct.rxnList{organismModelStruct.nRxns} = rxnNameWithCompartment;
organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn = organismModelStruct.nRxns; % at some point check that these are all used
organismModelStruct.rxns.(rxnNameWithCompartment).reactants = reactants; 
organismModelStruct.rxns.(rxnNameWithCompartment).products = products; 
organismModelStruct.rxns.(rxnNameWithCompartment).reversible = 0; 
organismModelStruct.rxns.(rxnNameWithCompartment).orrigString = 'see Biomass.txt file';
organismModelStruct.rxns.(rxnNameWithCompartment).rawRegulationString = '';
organismModelStruct.rxns.(rxnNameWithCompartment).metabolicClass = 'Biomass';
organismModelStruct.rxns.(rxnNameWithCompartment).metabolicSubClass = 'Biomass';
organismModelStruct.rxns.(rxnNameWithCompartment).rawGeneDependanceString = '';
% store reversibility
organismModelStruct.reversibility(organismModelStruct.nRxns) = 0;  
% store compartment
organismModelStruct.rxnCompartmentList{organismModelStruct.nRxns} = 'b'; % what copartment designation should biomass have?? (see a few lines above as well)

end % function parseBiomass

function organismModelStruct = addExchangeReactions(organismModelStruct)
% Add reactions to exchange 'e' () metabolites

% find 'e' metabolites
idxExchangeMetabolites = find( strcmp('e', organismModelStruct.metaboliteCompartmentList));

% loop thorugh exchange metabolites adding an exchange reaction
for iExchMet = 1 : length(idxExchangeMetabolites)
    
    % reaction name w/ exchange designation
    rxnNameWithCompartment = [organismModelStruct.metaboliteList{idxExchangeMetabolites(iExchMet)}];
    
    % update reaction counter
    organismModelStruct.nRxns = organismModelStruct.nRxns + 1;
    
    % save in structure
    organismModelStruct.rxnList{organismModelStruct.nRxns} = rxnNameWithCompartment;
    organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn = organismModelStruct.nRxns; % at some point check that these are all used
    organismModelStruct.rxns.(rxnNameWithCompartment).reactants = ...
        {organismModelStruct.metaboliteList{idxExchangeMetabolites(iExchMet)} 'e' [-1]}; 
    organismModelStruct.rxns.(rxnNameWithCompartment).products = {}; 
    organismModelStruct.rxns.(rxnNameWithCompartment).reversible = 1; 
    organismModelStruct.rxns.(rxnNameWithCompartment).orrigString = 'auto added exchange reaction';
    organismModelStruct.rxns.(rxnNameWithCompartment).rawRegulationString = '';
    organismModelStruct.mediaExchange.(rxnNameWithCompartment).idxMedia = iExchMet;
    organismModelStruct.mediaExchangeRxnIdxList(iExchMet) = ...
        organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn;
    organismModelStruct.rxns.(rxnNameWithCompartment).metabolicClass = 'Exchange(Model)';
    organismModelStruct.rxns.(rxnNameWithCompartment).metabolicSubClass = 'Exchange(Model)';
    organismModelStruct.rxns.(rxnNameWithCompartment).rawGeneDependanceString = '';

    % store reversibility
    organismModelStruct.reversibility(organismModelStruct.nRxns) = 1;  
    % store compartment
    organismModelStruct.rxnCompartmentList{organismModelStruct.nRxns} = 'e';
    
    % add to stoich matrix
    % add biomass reaction to stoich matrix
    rowStoichMatrix = organismModelStruct.metabolites.(organismModelStruct.metaboliteList{idxExchangeMetabolites(iExchMet)}).idxMetabolite;
    organismModelStruct.stoichMatrix(rowStoichMatrix, organismModelStruct.nRxns) = -1;
    
end % for iExchMet

organismModelStruct.nMediaExchange = length(idxExchangeMetabolites);

end % function addExchangeReactions

function organismModelStruct = addAccumulationReactions(organismModelStruct)
% Add reactions to accumulate 'a' metabolites

% find 'a' metabolites
idxNonSSMetabolites = find( strcmp('a', organismModelStruct.metaboliteCompartmentList));

% loop thorugh exchange metabolites adding an exchange reaction
for iNonSSMet = idxNonSSMetabolites
    
    % reaction name w/ exchange designation
    rxnNameWithCompartment = [organismModelStruct.metaboliteList{iNonSSMet}];
    
    % update reaction counter
    organismModelStruct.nRxns = organismModelStruct.nRxns + 1;
    
    % save in structure
    organismModelStruct.rxnList{organismModelStruct.nRxns} = rxnNameWithCompartment;
    organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn = organismModelStruct.nRxns; % at some point check that these are all used
    organismModelStruct.rxns.(rxnNameWithCompartment).reactants = {organismModelStruct.metaboliteList{iNonSSMet} 'a' [-1]}; 
    organismModelStruct.rxns.(rxnNameWithCompartment).products = {}; 
    organismModelStruct.rxns.(rxnNameWithCompartment).reversible = 0; 
    organismModelStruct.rxns.(rxnNameWithCompartment).orrigString = 'auto added Non SS reaction';
    organismModelStruct.rxns.(rxnNameWithCompartment).rawRegulationString = '';
    organismModelStruct.rxns.(rxnNameWithCompartment).metabolicClass = 'Accumulation';
    organismModelStruct.rxns.(rxnNameWithCompartment).metabolicSubClass = 'Accumulation';
    organismModelStruct.rxns.(rxnNameWithCompartment).rawGeneDependanceString = '';

    % store reversibility
    organismModelStruct.reversibility(organismModelStruct.nRxns) = 0;  
    % store compartment
    organismModelStruct.rxnCompartmentList{organismModelStruct.nRxns} = 'a';
    
    % add to stoich matrix
    rowStoichMatrix = organismModelStruct.metabolites.(organismModelStruct.metaboliteList{iNonSSMet}).idxMetabolite;
    organismModelStruct.stoichMatrix(rowStoichMatrix, organismModelStruct.nRxns) = -1;
    
end % for iNonSSMet

end % function addAccumulationReactions

function [organismModelStruct, organismModelErrors] = checkMetaboliteOrphans(organismModelStruct)
% check for metabolites that can either only be created or destroyed, and store reactions
% participated in by each reactant in model structure for easy reference later
% also stores in organismModel the reactions associated with each metabolite

organismModelErrors = {};

% for all metabolites
for iMetCheck = 1 : organismModelStruct.nMetabolites
    
    nameMetToCheck = organismModelStruct.metaboliteList{iMetCheck};
    
    % metabolites that can either only be created or only removed
    idxMetaboliteRow = organismModelStruct.metabolites.(nameMetToCheck).idxMetabolite;
    metaboliteRowOfStoichMatrix = organismModelStruct.stoichMatrix( idxMetaboliteRow ,:);
    
    % test for all positive or all negative stoichiometric coeffieicents
    allStoichPositive = all(metaboliteRowOfStoichMatrix >= 0);
    allStoichNegative = all(metaboliteRowOfStoichMatrix <= 0);
    
    % test for any reversibility in reactions containing met
    noReversibility = ~any( metaboliteRowOfStoichMatrix .* organismModelStruct.reversibility );
    
%     % add error if all of one sign and none are reversible
%     if (allStoichPositive || allStoichNegative) && noReversibility
%         organismModelErrors = [ organismModelErrors ;
%                                ['it was found that metabolite ', nameMetToCheck,...
%             ' is either only created or only removed from default model setup' ] ] ;
%     end % if all stoich of one sign
    
    % simulationously add a list of reactions each metabolite participates in to the
    % metabolite storage structure
    idxRxnsContainingMetabolite = ~(metaboliteRowOfStoichMatrix == 0);
    nameRxnsContainingMetabolite = organismModelStruct.rxnList{ idxRxnsContainingMetabolite };
    % store
    organismModelStruct.metabolites.(nameMetToCheck).rxnsUsingMet = nameRxnsContainingMetabolite;
    
end % for iMetCheck

% store all unique compartments
organismModelStruct.compartmentList = unique(organismModelStruct.metaboliteCompartmentList);

end % function checkMetaboliteOrphans 

function [organismModelStruct, organismModelErrors] = ...
    parseFluxBounds(rawTextReadFluxBounds, organismModelStruct, fluxBoundsDefault)

% default upper and lower bounds to use
fluxBoundDefaultUpper = fluxBoundsDefault(2);
fluxBoundDefaultLower = fluxBoundsDefault(1);

% before parsing for specific limits, use defaults and reversibility
organismModelStruct.fluxBoundsDefault = [organismModelStruct.reversibility' * fluxBoundDefaultLower , ...
                                         ones(organismModelStruct.nRxns,1) * fluxBoundDefaultUpper ] ; 
                                     

% separate into cells
rawTextReadFluxBoundsName = rawTextReadFluxBounds{ 1 };
rawTextReadFluxBoundsCompartment = rawTextReadFluxBounds{ 2 };
rawTextReadFluxBoundsLower = rawTextReadFluxBounds{ 3 };
rawTextReadFluxBoundsUpper = rawTextReadFluxBounds{ 4 };
rawTextReadFluxBoundsCondition = rawTextReadFluxBounds{ 5 };
if length( rawTextReadFluxBoundsCondition ) < length( rawTextReadFluxBoundsUpper  )
    rawTextReadFluxBoundsCondition{ length( rawTextReadFluxBoundsUpper  ) } = '';
end

% number of flux bounds read
nFluxBounds = size( rawTextReadFluxBoundsName, 1);

% intialize error storage
organismModelErrors = {};

nFluxBoundsConditional = 0; % counter for conditional flux bounds
organismModelStruct.fluxBoundsConditional = {};


% loop through setting flux bounds
for iFluxBound = 1 : nFluxBounds
    
    rxnNameWithCompartment = [lower(rawTextReadFluxBoundsCompartment{iFluxBound})...
                              upper(rawTextReadFluxBoundsName{iFluxBound})];
    
    % check that limited flux is represented in model
    if isfield(organismModelStruct.rxns, rxnNameWithCompartment)
        
       
        nonWhitespaceCondition = regexp(  rawTextReadFluxBoundsCondition{iFluxBound}, ' \S', 'match');
        
        if isempty( nonWhitespaceCondition ) % || ~ organismModelStruct.regulationFlag
        
            % flux found in model, set upper and lower bound
            idxFluxToBound = organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn;
            fluxBoundLower = rawTextReadFluxBoundsLower(iFluxBound);
            fluxBoundUpper = rawTextReadFluxBoundsUpper(iFluxBound);

            % verify that flux upper bound is equal to or greater than lower limit
            if fluxBoundUpper >= fluxBoundLower

                % upper bound is greater, store
                organismModelStruct.fluxBoundsDefault(idxFluxToBound,:) = [fluxBoundLower fluxBoundUpper];

            else

                % bounds inconsistant, error
                organismModelErrors = [ organismModelErrors ;
                                     ['flux upper bound is less than lower bound for: ',...
                                     rxnNameWithCompartment] ]; 

            end % if upper/lower bounds consistant
            
        else
            % non whitespace found for condition
            
            % flux found in model, set upper and lower bound
            idxFluxToBound = organismModelStruct.rxns.(rxnNameWithCompartment).idxRxn;
            fluxBoundLower = rawTextReadFluxBoundsLower(iFluxBound);
            fluxBoundUpper = rawTextReadFluxBoundsUpper(iFluxBound);

            % verify that flux upper bound is equal to or greater than lower limit
            if ~ (fluxBoundUpper >= fluxBoundLower)

                % bounds inconsistant, error
                organismModelErrors = [ organismModelErrors ;
                                     ['flux upper bound is less than lower bound for: ',...
                                     rxnNameWithCompartment] ]; 
            else
                
                [boundRegulationRulesEvalString, ~, boundRegulationErrors ] = ...
                FBA_EB_regulationRuleParse( rawTextReadFluxBoundsCondition{iFluxBound} , rxnNameWithCompartment, organismModelStruct  );
                
                if isempty( boundRegulationErrors )
                    
                    nFluxBoundsConditional = nFluxBoundsConditional + 1;
                    
                    organismModelStruct.fluxBoundsConditional{ nFluxBoundsConditional, 1 } = idxFluxToBound;
                    organismModelStruct.fluxBoundsConditional{ nFluxBoundsConditional, 2 } = [fluxBoundLower, fluxBoundUpper];
                    organismModelStruct.fluxBoundsConditional{ nFluxBoundsConditional, 3 } = boundRegulationRulesEvalString;
                    
                else
                    
                    organismModelErrors = [ organismModelErrors ; boundRegulationErrors];
                    
                end
                              

            end % if upper/lower bounds inconsistant
        
        end % if no condition
        
    else
        % flux whose limit is being set does not exist in model
        organismModelErrors = [ organismModelErrors ;
                                 ['flux whos bound is being set: ',...
                                 rxnNameWithCompartment,' is not found in the model'] ]; 
    end % if flux (bounded) is found in model
    
end % for iFluxBounds

end % function parseFluxBounds

function [organismModelStruct, organismModelErrors] = parseGenes(rawTextReadGenes, organismModelStruct, geneAutoAddOptions)

% separate into columns
rawTextReadGenesNumber = rawTextReadGenes{ 1 };
rawTextReadGenesName = rawTextReadGenes{ 2 };
rawTextReadGenesRule = rawTextReadGenes{ 3 };
rawTextReadGenesDelay = rawTextReadGenes{ 4 };

% numer of genes determined by length of imported data
nGene = length( rawTextReadGenesNumber );

% initialize storage of specified vs full whole numbering of e coli genome
allGene = zeros( geneAutoAddOptions.maxGeneNumber, 1); 

% initialize error storage
organismModelErrors = {};

% for all genes in raw text add to storage structure with raw rule, and idx
for iGene = 1 : nGene
    
    % get gene number
    rawGeneNumber = rawTextReadGenesNumber{ iGene };
    
    % check that gene number matches basic format
    [~, ~, ~, geneNumMatch, ~, ~, geneNumSplit] = regexp( rawGeneNumber, 'b\d{4}', 'once' );
    [~, ~, ~, unexpectedMatch, ~, ~, ~] = regexp( cat(2, geneNumSplit{:}), '\S', 'once' );
    
    % if match found and no unexpected (other than whitespace) extras, store
    if ~isempty(geneNumMatch) && isempty(unexpectedMatch)
        
        geneNumber = geneNumMatch;
        geneNumberOnly = str2double( geneNumber(2:end) );
        organismModelStruct.genesList{iGene} = geneNumber;
        organismModelStruct.genes.(geneNumber).geneName = rawTextReadGenesName{ iGene };
        organismModelStruct.genes.(geneNumber).rawRegulationString = rawTextReadGenesRule{ iGene };
        organismModelStruct.genes.(geneNumber).idxGene = iGene ;
        organismModelStruct.genes.(geneNumber).geneNumberOnly = geneNumberOnly;
        organismModelStruct.genes.(geneNumber).activityDelay = rawTextReadGenesDelay(iGene);
        organismModelStruct.activityDelayList(iGene) = rawTextReadGenesDelay(iGene);
        allGene(geneNumberOnly) = 1;
        
    else
        organismModelErrors = [ organismModelErrors ;
                                 {['gene number string: ''',rawGeneNumber,...
                                 ''' is improperly formatted or contains unexpected characters']} ]; 
    
    end % if only expected gene number found in string
    
end % for iGene

% if desired, auto add all whole numbers not user specified
geneAutoAddFlag = geneAutoAddOptions.flag;  
activityDelayDefault = geneAutoAddOptions.delayDefault; 

if geneAutoAddFlag
    
    geneNumbersAutoAdd = find( ~allGene );
    nGeneAutoAdd = length( geneNumbersAutoAdd );
    
    for iGeneAutoAdd = 1 : nGeneAutoAdd
        
        % create pad of zeros for blattner number auto gen
        if geneNumbersAutoAdd(iGeneAutoAdd) > 999
            strPad = '';
        elseif geneNumbersAutoAdd(iGeneAutoAdd) > 99
            strPad = '0';
        elseif geneNumbersAutoAdd(iGeneAutoAdd) > 9
            strPad = '00';
        else
            strPad = '000';
        end % if switch for string padding
        
        geneNumber = ['b',strPad,num2str(geneNumbersAutoAdd(iGeneAutoAdd))];
        geneNumberOnly = geneNumbersAutoAdd(iGeneAutoAdd);
        organismModelStruct.genesList{ nGene + iGeneAutoAdd  } = geneNumber;
        organismModelStruct.genes.(geneNumber).geneName = [geneNumber,'autoAdded'];
        organismModelStruct.genes.(geneNumber).rawRegulationString = '';
        organismModelStruct.genes.(geneNumber).idxGene = nGene + iGeneAutoAdd ;
        organismModelStruct.genes.(geneNumber).geneNumberOnly = geneNumberOnly;
        organismModelStruct.genes.(geneNumber).activityDelay = activityDelayDefault;
        organismModelStruct.activityDelayList(nGene + iGeneAutoAdd) = activityDelayDefault;
        
    end % for iGeneAutoAdd
    
    % update total gene count
    nGene = nGene + nGeneAutoAdd;
    
end % if geneAutoAddFlag 

% store final count of genes
organismModelStruct.nGene = nGene;

end % function parseGenes

function [organismModelStruct, organismModelErrors] = parseTFs(rawTextReadTFs, organismModelStruct)

% initialize error storage
organismModelErrors = {};

% separate into columns
rawTextReadTFsName = rawTextReadTFs{ 1 };
rawTextReadTFsRule = rawTextReadTFs{ 2 };

% catch for uneven column lengths if no rule for last entry
if size(rawTextReadTFsName,1) - size(rawTextReadTFsRule,1) == 1;
    % one less entry found for reaction rules than reaction name/eqns, assume this was an
    % empty cell not captured in tab delimiting of text and add a whitespace string there
    rawTextReadTFsRule = [rawTextReadTFsRule; ' '];
end % if one less row in regulation (blank last cell)

% numer of TF determined by length of imported data
nTFs = length( rawTextReadTFsName );

% for all TF in raw text add to storage structure with raw rule, and idx
for iTF = 1 : nTFs
    
    % get raw TF name
    rawTFname = rawTextReadTFsName{ iTF };
    
    % check that TF name is valid variable
    [~, ~, ~, TFnameMatch, ~, ~, TFnameSplit] = regexp( rawTFname, '[a-zA-Z]\w+', 'once' );
    reservedKeywordFound = FBA_EB_reservedKeywordCheck(rawTFname);
    [~, ~, ~, unexpectedMatch, ~, ~, ~] = regexp( cat(2, TFnameSplit{:}), '\S', 'once' );
    
    % if match found and no unexpected (other than whitespace) extras, store
    if ~isempty(TFnameMatch) && isempty(unexpectedMatch) && ~reservedKeywordFound
        
        organismModelStruct.tfList{iTF} = TFnameMatch;
        organismModelStruct.TFs.(TFnameMatch).geneName = TFnameMatch;
        organismModelStruct.TFs.(TFnameMatch).rawRegulationString = rawTextReadTFsRule { iTF };
        organismModelStruct.TFs.(TFnameMatch).idxTF = iTF ;
        
    else
        organismModelErrors = [ organismModelErrors ;
                                 {['TF name string: ''',TFnameMatch,...
                                 ''' is improperly formatted or contains unexpected characters']} ]; 
    
    end % if only expected TF number found in string
    
end % for iTF

% store number of TFs
organismModelStruct.nTFs = nTFs;

end % function parseTFs

function [organismModelStruct, stimuliErrors] = parseStimuli(rawTextReadStimuli, organismModelStruct)

stimuliErrors = {};

rawTextReadStimuliNames = rawTextReadStimuli{ 1 };
rawTextReadStimuliValues = rawTextReadStimuli{ 2 };

% dummy field name, removed later
organismModelStruct.stimuli.dummy.idxStim = 0;

% number of stimuli found
organismModelStruct.nStimuli = size(rawTextReadStimuliNames, 1);

% loop through stimuli storing in structure
for iStim = 1 : organismModelStruct.nStimuli
    
    % check validity of name, after making caps and replacing any dashes with underscores
    stimName = upper(strrep(rawTextReadStimuliNames{ iStim },'-','_'));
    
    if isvarname(stimName) && ~FBA_EB_reservedKeywordCheck(stimName) ...
            && ~isfield(organismModelStruct.stimuli, stimName )
        
        % name is valid matlab variable, does not contain reserved keywords, and is not
        % already present in stimuli list
        % store in list, etc
        organismModelStruct.stimuliList{iStim} = stimName;
        organismModelStruct.stimuli.(stimName).idxStim = iStim;
        
        % store value as either logical or
        stimuliValueLogical = {'on', [true];
                               'true', [true];
                               'false', [false];
                               'off', [false]};
        stimuliValueLogicalTest = strcmpi( rawTextReadStimuliValues{iStim},  stimuliValueLogical(:,1) );
        if any(stimuliValueLogicalTest)
            
            % logical on/off string is found, replace with true/false string
            organismModelStruct.stimuliValuesDefault{iStim} = stimuliValueLogical{ stimuliValueLogicalTest, 2 };
            organismModelStruct.stimuli.(stimName).valueDefault = stimuliValueLogical{ stimuliValueLogicalTest, 2 };
            
        elseif ~isnan( str2double( rawTextReadStimuliValues{iStim} ) )
            
            % number is found
            organismModelStruct.stimuliValuesDefault{iStim} = str2double( rawTextReadStimuliValues{iStim} );
            organismModelStruct.stimuli.(stimName).valueDefault = str2double( rawTextReadStimuliValues{iStim} );
            
        else
            
            % value is niether logical or numeric
            stimErrors = [ stimErrors ; ...
               {['stimuli value: ', rawTextReadStimuliValues{iStim}, 'associated with ',stimName,', is not recognized as logical or numeric' ]}];
                    
        end % if for stim values
        
            
    else
        % name is NOT valid matlab variable, DOES contain reserved keywords, or IS
        % already present in stimuli list
        stimErrors = [ stimErrors ; ...
           {['stimuli name: ', stimName,', is invalid, contains a keyword, or has already been used' ]}];
    end

end %

end % function parseStimuli

function [organismModelStruct, organismModelErrors] = parseTFsRules(organismModelStruct)

organismModelErrors = {};

organismModelStruct.tfRegulationRulesEvalStrings{organismModelStruct.nTFs} = '';

for iTFsReg = 1 : organismModelStruct.nTFs
    
    % name of reaction being considered this iteration
    tfName = organismModelStruct.tfList{ iTFsReg };
    tfRawRegulationString = organismModelStruct.TFs.(tfName).rawRegulationString;
    
    % parse rules using function FBA_EB_rxnRegulationRuleParse.m
    [tfRegulationRulesEvalString, tfDependancies, tfRegulationErrors ] = ...
                FBA_EB_regulationRuleParse( tfRawRegulationString, tfName, organismModelStruct );
    % store gene dependances in rxn
    organismModelStruct.tfRegulationRulesEvalStrings{iTFsReg} = tfRegulationRulesEvalString;
    organismModelStruct.TFs.(tfName).tfDependancies = tfDependancies;
    organismModelStruct.TFs.(tfName).tfRegulationRulesEvalString = tfRegulationRulesEvalString;
    
    %
    organismModelErrors = [ organismModelErrors ; tfRegulationErrors];
    
end % for iTFsReg

end % function parseTFsRules

function [organismModelStruct, organismModelErrors] = parseGenesRules(organismModelStruct)

% initialize storage variables
organismModelErrors = {};
organismModelStruct.geneRegulationRulesEvalStrings{organismModelStruct.nGene} = '';

for iGeneReg = 1 : organismModelStruct.nGene
    
    % name of reaction being considered this iteration
    geneName = organismModelStruct.genesList{ iGeneReg };
    geneRawRegulationString = organismModelStruct.genes.(geneName).rawRegulationString;
    
    % parse rules using function FBA_EB_rxnRegulationRuleParse.m
    [geneRegulationRulesEvalString, geneDependancies, geneRegulationErrors ] = ...
                FBA_EB_regulationRuleParse( geneRawRegulationString, geneName, organismModelStruct );
    % store gene dependances in rxn
    organismModelStruct.geneRegulationRulesEvalStrings{iGeneReg} = geneRegulationRulesEvalString;
    organismModelStruct.genes.(geneName).geneDependancies = geneDependancies;
    organismModelStruct.genes.(geneName).geneRegulationRulesEvalString = geneRegulationRulesEvalString;
    
    %
    organismModelErrors = [ organismModelErrors ; geneRegulationErrors];
    
end % for iGenesReg

end % function parseGenesRules

function [organismModelStruct, organismModelErrors] = parseRxnsRules(organismModelStruct)

organismModelErrors = {};

organismModelStruct.rxnRegulationRulesEvalStrings{organismModelStruct.nRxns} = '';

for iRxnReg = 1 : organismModelStruct.nRxns
    
    % name of reaction being considered this iteration
    rxnName = organismModelStruct.rxnList{ iRxnReg };
    rxnRawRegulationString = organismModelStruct.rxns.(rxnName).rawRegulationString;
    % parse rules using function FBA_EB_rxnRegulationRuleParse.m
    [rxnRegulationRulesEvalString, rxnGeneDependancies, rxnRegulationErrors ] = ...
                FBA_EB_regulationRuleParse( rxnRawRegulationString, rxnName, organismModelStruct  );
%     [rxnRegulationRulesEvalString, rxnGeneDependancies, rxnRegulationErrors ] = FBA_EB_rxnRegulationRuleParse( rxnName, organismModelStruct );
    
    % store gene dependances in rxn
    organismModelStruct.rxnRegulationRulesEvalStrings{iRxnReg} = rxnRegulationRulesEvalString;
%     organismModelStruct.rxns.(rxnName).rxnGeneDependancies = rxnGeneDependancies;
    organismModelStruct.rxns.(rxnName).rxnRegulationRulesEvalString = rxnRegulationRulesEvalString;
    
    %
    organismModelErrors = [ organismModelErrors ; rxnRegulationErrors];
    
end % for iRxnReg

end % function parseRxnsRules

function [organismModelStruct, organismModelErrors] = parseRxnsGenes(organismModelStruct)

organismModelErrors = {};

organismModelStruct.rxnGeneDependance{organismModelStruct.nRxns} = '';

for iRxnReg = 1 : organismModelStruct.nRxns
    
    % name of reaction being considered this iteration
    rxnName = organismModelStruct.rxnList{ iRxnReg };
    rawGeneDependanceString = organismModelStruct.rxns.(rxnName).rawGeneDependanceString;
    % parse rules using function FBA_EB_rxnRegulationRuleParse.m
    [rxnGeneDependanceEvalStrings, ~, rxnRegulationErrors ] = ...
                FBA_EB_regulationRuleParse( rawGeneDependanceString, rxnName, organismModelStruct  );
%     [rxnRegulationRulesEvalString, rxnGeneDependancies, rxnRegulationErrors ] = FBA_EB_rxnRegulationRuleParse( rxnName, organismModelStruct );
    
    % store gene dependances in rxn
    organismModelStruct.rxnGeneDependanceEvalStrings{iRxnReg} = rxnGeneDependanceEvalStrings;
    organismModelStruct.rxns.(rxnName).rxnGeneDependanceEvalStrings = rxnGeneDependanceEvalStrings;
    
    %
    organismModelErrors = [ organismModelErrors ; rxnRegulationErrors];
    
end % for iRxnReg

end % function parseRxnsRules

function [organismModelStruct, unnamedMet] = parseMetaboliteNames(rawTextReadMetabolites, organismModelStruct)

organismModelWarnings = {};

rawTextReadMetabolitesAbrev = rawTextReadMetabolites{ 1 };
rawTextReadMetabolites  = rawTextReadMetabolites{ 2 };

nMetName = length(rawTextReadMetabolitesAbrev);
nCompartments = length( organismModelStruct.compartmentList );

% initialize storage of names of metabolites ( same as names so that any met w/o
% abbreviations are 'named' their abreviation )
organismModelStruct.metaboliteNames = organismModelStruct.metaboliteList;
logicalNammedMetabolites = false( organismModelStruct.nMetabolites, 1 );

for iMetName = 1 : nMetName
    
    % convert abreviation to form used in model struct (cap, dash to underscore)
    abrevConverted = upper( strrep( rawTextReadMetabolitesAbrev{1}, '-', '_' ) );
    
    % for each compartment find matching metabolite abreviations
    for iComp = 1 : nCompartments
        
        logicalMet = strcmp( [ organismModelStruct.compartmentList{ iComp }, abrevConverted ],...
                        organismModelStruct.metaboliteList );
                    
        if any(logicalMet)
            
            % store full metabolite name
            organismModelStruct.metaboliteNames{ logicalMet } = rawTextReadMetabolites{ iMetName };
            
            % flag that metabolite as having name
            logicalNammedMetabolites( logicalMet ) = true;
            
        end % if any logicalMet
        
    end % for iComp
    
end % for iMetName

unnamedMet = organismModelStruct.metaboliteList( ~logicalNammedMetabolites )';

end % function parseMetaboliteNames

function [organismModelStruct ] = parseMetabolicClass( organismModelStruct )

% create storage
metabolicClassesRxns = repmat({''},organismModelStruct.nRxns,1); 
metabolicClassesIdx = zeros(organismModelStruct.nRxns,1);  

for iMetClass = 1 : organismModelStruct.nRxns
    
    metabolicClass = organismModelStruct.rxns.(organismModelStruct.rxnList{ iMetClass }).metabolicClass;
    % find any
    firstExistingMatch = find( strcmpi(metabolicClass, metabolicClassesRxns), 1, 'first');
    if ~isempty( firstExistingMatch )
        % replace with existing case insensitiv e match to apply capitilzation homegenaity
        organismModelStruct.rxns.(organismModelStruct.rxnList{ iMetClass }).metabolicClass = ...
            metabolicClassesRxns{ firstExistingMatch };
        metabolicClass = metabolicClassesRxns{ firstExistingMatch };
    end % end if existing match
    
    % store in list of all rxns metabolic classes
    metabolicClassesRxns{ iMetClass } = metabolicClass;
    
end % for iMetClass

% create unique list and idx list referencing this
metabolicClassList = unique( metabolicClassesRxns );
nMetabolicClass = length( metabolicClassList );

for iUnMetClass = 1 : nMetabolicClass
    
    % find all rxns with this class
    logicalRxnsThisClass = strcmp( metabolicClassList{ iUnMetClass }, metabolicClassesRxns );
    
    % store idx ref to this class for each rxn
    metabolicClassesIdx( logicalRxnsThisClass ) = repmat( iUnMetClass, sum( logicalRxnsThisClass ), 1);
    
end % for iUnMetClass 

% store it all

organismModelStruct.metabolicClassList = metabolicClassList;
organismModelStruct.nMetabolicClass = nMetabolicClass;
organismModelStruct.metabolicClassesRxns = metabolicClassesRxns; 
organismModelStruct.metabolicClassesIdx = metabolicClassesIdx ; 


% now parse subclasses
% create storage
metabolicSubClassesRxns = repmat({''},organismModelStruct.nRxns,1); 

for iMetSubClass = 1 : organismModelStruct.nRxns
    
    metabolicSubClass = organismModelStruct.rxns.(organismModelStruct.rxnList{ iMetSubClass }).metabolicSubClass;
    % find any
    firstExistingMatch = find( strcmpi(metabolicSubClass, metabolicSubClassesRxns), 1, 'first');
    if ~isempty( firstExistingMatch )
        % replace with existing case insensitiv e match to apply capitilzation homegenaity
        organismModelStruct.rxns.(organismModelStruct.rxnList{ iMetSubClass }).metabolicSubClass = ...
            metabolicSubClassesRxns{ firstExistingMatch };
        metabolicSubClass = metabolicSubClassesRxns{ firstExistingMatch };
    end % end if existing match
    
    % store in list of all rxns metabolic classes
    metabolicSubClassesRxns{ iMetSubClass } = metabolicSubClass;
    
end % for iMetSubClass

% loop through classes identifying unique subclasses within
metabolicSubClasses{nMetabolicClass,2} = [];
for iClass = 1 : nMetabolicClass
    
    reactionSubClasses = metabolicSubClassesRxns( metabolicClassesIdx == iClass );
    uniqueSubClasses = unique( reactionSubClasses );
    metabolicSubClasses{ iClass, 1} = uniqueSubClasses;
    
    for iSubClass = 1 : length( uniqueSubClasses );
        
        metabolicSubClasses{ iClass, 2}{ iSubClass } = ...
            find( strcmp( uniqueSubClasses{iSubClass},reactionSubClasses) );
        
    end % for iSubClass
    
end % for iClass

organismModelStruct.metabolicSubClasses = metabolicSubClasses;

end % function parseMetabolicClass

function [ organismModelStruct ] = findSimpleCycles(organismModelStruct)

% loop through all of organism rxns determining those that are equivalent
% stoichiometrically and therefore likely to be cycles durring simulation

organismModelStruct.potentialCycles = {};
nCycles = 0;
alreadyMatched = false( organismModelStruct.nRxns, 1 );

for iRxns = 1 : organismModelStruct.nRxns
    if ~alreadyMatched( iRxns )
        idxCycles = [iRxns];
        flagCycle = false;
        for jRxns = iRxns + 1 : organismModelStruct.nRxns
            % test for stoichiometric eqivalency
            sameForward = organismModelStruct.stoichMatrix(:,iRxns) == organismModelStruct.stoichMatrix(:,jRxns);
            sameReverse = organismModelStruct.stoichMatrix(:,iRxns) == - organismModelStruct.stoichMatrix(:,jRxns);
            % store under various conditions of cycle potential
            % same forward means at least one must be reversible for cycle to occur
            if all(sameReverse) || (all(sameForward) && ( organismModelStruct.reversibility(iRxns) || ...
                    organismModelStruct.reversibility(jRxns) ))

                idxCycles = [ idxCycles, jRxns ];

                flagCycle = true;
                
                alreadyMatched( jRxns ) = true;

            end % if rxn j is a cycles with rxn i
        end % for jRxns

        if flagCycle
            nCycles = nCycles + 1;
            organismModelStruct.potentialCycles{nCycles,1} = idxCycles;
        end % if cycle exists to store
    end
end % for iRxns

end % findSimpleCycles

function [ organismModelStruct ] = parseCycles( organismModelStruct ,rawTextReadCycles )


    nCycRaw = length( rawTextReadCycles{1});
    rxnListStringsRaw = rawTextReadCycles{1};
    
    nCyc = 0;
    
    idxRxnsCycle = [];
    
    for iCycRaw = 1:nCycRaw
        
        rxnString = rxnListStringsRaw{ iCycRaw };
            
            rxnString = strrep(rxnString,' ','');
            rxnString = strrep(rxnString,'"','');
            rxnString = strrep(rxnString,'-','_');
            rxnString = [lower(rxnString(1)), upper(rxnString(2:end))];
            
            if isfield( organismModelStruct.rxns, rxnString)
                nCyc = nCyc + 1;
                idxRxnsCycle( nCyc ) = organismModelStruct.rxns.( rxnString ).idxRxn;
            else
                disp(['rxn: ',rxnString,' , listed in cycles not found.'])
            end
        
    end % for iCycRaw
    
    organismModelStruct.cycleIdxsUser = idxRxnsCycle;
%     organismModelStruct.cycleCoeffsUser = cycleCoeffs;
%     

%     nCycRaw = length( rawTextReadCycles{1});
%     rxnListStringsRaw = rawTextReadCycles{1};
%     coeffListStringsRaw = rawTextReadCycles{2};
%     
%     nCyc = 0;
%     
%     cycleIdxs = {};
%     cycleCoeffs = {};
%     
%     for iCycRaw = 1:nCycRaw
%         
%         rxnListString = rxnListStringsRaw{ iCycRaw }
%         coeffListString = coeffListStringsRaw{ iCycRaw }
%         
%         idxCommasRxnList = strfind( rxnListString, ',' );
%         idxCommasCoeffList = strfind( coeffListString, ',' );
%         
%         nCycRxns = length( idxCommasRxnList )+1;
%         
%         idxRxnsCycle = zeros(1,nCycRxns);
%         coeffRxnsCycle = zeros(1,nCycRxns);
%         
%         for iCycRxns = 1:nCycRxns
%            
%             if iCycRxns == 1
%                 rxnString = rxnListString(1:idxCommasRxnList(iCycRxns)-1);
%                 coeffString = coeffListString(1:idxCommasCoeffList(iCycRxns)-1);
%             elseif iCycRxns == nCycRxns
%                 rxnString = rxnListString(idxCommasRxnList(iCycRxns-1)+1:end);
%                 coeffString = coeffListString(idxCommasCoeffList(iCycRxns-1)+1:end);
%             else
%                 rxnString = rxnListString(idxCommasRxnList(iCycRxns-1)+1:idxCommasRxnList(iCycRxns)-1);
%                 coeffString = coeffListString(idxCommasCoeffList(iCycRxns-1)+1:idxCommasCoeffList(iCycRxns)-1);
%             end % if first of last
%             
%             rxnString = strrep(rxnString,' ','');
%             rxnString = strrep(rxnString,'"','');
%             rxnString = strrep(rxnString,'-','_');
%             coeffString = strrep(coeffString,' ','');
%             coeffString = strrep(coeffString,'"','');
%             
%             rxnString = [lower(rxnString(1)), upper(rxnString(2:end))];
%             
%             coeffRxnsCycle(iCycRxns) = str2double( coeffString );
%             
%             if isfield( organismModelStruct.rxns, rxnString)
%                 idxRxnsCycle( iCycRxns ) = organismModelStruct.rxns.( rxnString ).idxRxn;
%             else
%                 disp(['rxn: ',rxnString,' , listedin cycles not found.'])
%             end
%             
%         end % for iCycRxns
%         
%         
%         % check that rxns indeed cycle
%         cycleRxnStoich = organismModelStruct.stoichMatrix(:,idxRxnsCycle);
%         coeffMatrix = diag( coeffRxnsCycle );
%         cycleRxnNet = sum( (cycleRxnStoich * coeffMatrix) , 2);
%         nonzeroNet = any(cycleRxnNet);
%         if nonzeroNet
%             disp(['nonzero net through cycle, in: ',num2str(iCycRaw)])
%         else
%             nCyc = nCyc + 1;
%             cycleIdxs{ nCyc } = idxRxnsCycle;
%             cycleCoeffs{ nCyc } = coeffRxnsCycle;
%         end
%         
%     end % for iCycRaw
%     
%     organismModelStruct.cycleIdxsUser = cycleIdxs;
%     organismModelStruct.cycleCoeffsUser = cycleCoeffs;
%     

end % parse cycles