% 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: After evaluation of FBA linear optimization at a given timestep, update the concentrations of biomass, 
% 				reset or consumed media, and deal wih any metabolite pooling.
% 
%
% Inputs: 	organism - organism structure with simulation fluxes detrmined for this siultation itteration 
%			iSim -  simulation itteration
% 			tStep - timestep length in hours
%
% Outputs: 	organismUpdated - organism structure iwth biomass concentration etc updated
% 			mediaExchanged - media consumed by this organism
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function [organismUpdated, mediaExchanged] = FBA_EB_updateGrowthAndExchange(organism, iSim, tStep)

% store growth           
organism.simulation.growth(iSim) = ...
    organism.simulation.fluxes(organism.model.rxns.aBIOMASS.idxRxn, iSim);


% store accumulation
organism.simulation.accumulation(:,iSim) = ...
    organism.model.stoichMatrix * organism.simulation.fluxes(:,iSim);

% update biomass concentration
organism.simulation.biomass( iSim ) = ...
    organism.simulation.biomass( iSim - 1) * exp(organism.simulation.growth(iSim) * tStep);

% calculate media exchanges, and metabolite pools
% test for growth, calc exchange (to avoid devide by zero)
if (organism.simulation.growth(iSim) > 0)

    %
    %
    %
    mediaExchanged = - organism.simulation.fluxes(organism.model.mediaExchangeRxnIdxList,iSim) ...
                                * organism.simulation.biomass(iSim-1) ...
                                * (1 - exp(tStep * organism.simulation.growth(iSim)))  ...
                                / organism.simulation.growth(iSim);
    
    %
    %
    %
    organism.simulation.metabolitePools(:, iSim) = organism.simulation.metabolitePools(:, iSim - 1) ...
        + organism.simulation.accumulation(:,iSim) * organism.simulation.biomass(iSim-1) ...
                                * (exp(tStep * organism.simulation.growth(iSim)) - 1)  ...
                                / organism.simulation.growth(iSim);
else
    % This is the case when there is no growth, jut uptake/out from current biomass
    %
    %
    mediaExchanged = organism.simulation.fluxes(organism.model.mediaExchangeRxnIdxList,iSim) ...
                                * organism.simulation.biomass(iSim-1) * tStep;
                            
    %
    %
    organism.simulation.metabolitePools(:, iSim) = organism.simulation.metabolitePools(:, iSim - 1) ...
        + organism.simulation.accumulation(:,iSim) * organism.simulation.biomass(iSim-1) * tStep;

end % if growth

% round any small negatives arising in metabolite pooling
organism.simulation.metabolitePools(:, iSim) = max( 0 , organism.simulation.metabolitePools(:, iSim));

% clean up cycles
organism.simulation.fluxes(:,iSim) = FBA_EB_cycleRemove(organism, iSim);

organismUpdated = organism ;