% 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: For one step of an uninfected organism calls functions to run FBA problem. 
% 			For infected organism evaluates rules and sets flux bounds, then runs virus and applies to the host.
% 
%
% Inputs: 	tStep - timestep in hours
%			iSim - current itteraction
%			organism - organism structure
%			mediaAvalible - media vector of avalible resources for this organism
%			stimuliCurrent - (placeholder for futer inplementation of stimuli inputs)
%
% Outputs:  organismSimulation - updated state of organism
%			mediaExchanged - media consumed by this organism this itteration 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: > stimuli is currently not implemented fully
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function [organismSimulation, mediaExchanged] = ...
    FBA_EB_modularSingleStepPremtiveLimit(tStep, iSim, organism, mediaAvalible, stimuliCurrent )


% set host bounds
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

fluxConstraintFunction = organism.methods.fluxConstraintHandle;
%
[organism.simulation] = fluxConstraintFunction(tStep, iSim, organism,...
                            mediaAvalible(organism.mediaIdxTranslation),stimuliCurrent);
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 


% set unconstrained viral bounds (if infected), and rhs
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

if organism.infectionFlag
    % organism is infected

    % store a copy of the organism before any changes are made, (this may be unnesceary, check)
    % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    organismRevert = organism;
    organism.model = organism.modelInfected;
    % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

    
    % no limits on viral ODEs this run    
    viralConstraintFunction = organism.methods.viralConstraintHandle;
    [ viralODElimits, viralIdealRequest, warningStr, hostFluxes ] = viralConstraintFunction(tStep, iSim, organismRevert, mediaAvalible, stimuliCurrent);
    organism.simulation.viralLimits(:,iSim) = viralODElimits{2}; % save for debugging later
    organism.simulation.viralIdeal(:,iSim) = viralIdealRequest; % save for debugging
    organism.simulation.warningOut = [ organism.simulation.warningOut ; warningStr ];
    organism.simulation.hostLim(:,iSim) = hostFluxes ;
    
    % run viral ODEs and use to constrain host
    [ organism.simulation ] = FBA_EB_virusRunAndApplyToHost( organism , viralODElimits, iSim, tStep );
    
    % steady State Right Hand Side, except viral metabolites
    % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    viralMetaboliteIndicies = (organism.model.nHostMetabolites + 1 : organism.model.nMetabolites)' ;
    freeRhsViral = [ viralMetaboliteIndicies, zeros(size(viralMetaboliteIndicies))];
    accumulationFunction = organism.methods.accumulationSetHandle;
    lowBoundRhs = accumulationFunction(organism, tStep, iSim);
    rhsSpecifications = {[], freeRhsViral, [] ,lowBoundRhs ,[]}; % set (not zero), free, upper, lower, or double bounds set
    % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    
    
else
    % organism is uninfected, only set rhs to steady state
    rhsSpecifications = {[], [], [] ,[] ,[]}; %  set (not zero), free, upper, lower, or double bounds set
end


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

% evaluate model
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% biomass max objective function
objectiveCoefficientVector = zeros( organism.model.nRxns, 1 );
objectiveCoefficientVector(  organism.model.cycleIdxsUser ) = ones(size(organism.model.cycleIdxsUser) )*(-1e-5);
objectiveCoefficientVector(organism.model.rxns.aBIOMASS.idxRxn) = 1;

% simulate using linear programming optimization (currently glpk)
[organism.simulation.fluxes(:,iSim) ,organism.simulation.status(iSim), organism.simulation.dual(:,iSim)] = ... %  ((fba)met / hour)
    FBA_EB_glpkSolve( organism.model.stoichMatrix,...
                organism.simulation.fluxBoundsLower(:,iSim),...
                organism.simulation.fluxBoundsUpper(:,iSim),...
                objectiveCoefficientVector,...
                rhsSpecifications);

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% call function to set growth rate etc from fluxes and calculate net media exchanged
[organism, mediaExchangedOrganism] = FBA_EB_updateGrowthAndExchange(organism, iSim, tStep);
% translate from organism order(indicies) to whole simulation order for
% media(metabolite) exchange
mediaExchanged = zeros( size(mediaAvalible) );
mediaExchanged(organism.mediaIdxTranslation) = mediaExchangedOrganism;

% assign output of only simulation state
organismSimulation = organism.simulation;
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    


