% Author: Elsa W. Birch, ewbirch@stanford.edu, Covert Lab - Stanford Bioengineering
% Date: July 27th 2012
% Version: associated with production of the figures for Birch et al. 2012
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Summary: Runs viral ODEs with limits applied, then converts dynamics to FBA units steady flux and applies as bounds to the FBA problem
% 
% Inputs: 	tStep - timestep in hours
%  			iSim - iteration of simulation
%			organism - organism structure that includes an infected organism
%			viralODElimits - vector of ODE rate limits
%
% Outputs: organismSimulation - updated simulation with viral demand applied as bounds to the host FBA problem
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function [ organismSimulation ] = FBA_EB_virusRunAndApplyToHost( organism , viralODElimits, iSim, tStep )

% prevent rounding ror issues
viralODElimits{2} = viralODElimits{2} * 0.999 ;

% evaluate viral ODEs 
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
% prepare function handle
viralODEfn = organism.modelInfected.handleViralODEs;
% prepare time variable range
tViralLim = ((iSim-2)*tStep: tStep/organism.simulation.nViralSubSteps :(iSim-1)*tStep) * 60 * 60; % hour to sec (FBA time to virus time)
% viral initial conditions
cViralInitial = organism.simulation.viralConcentrations( : , iSim - 1 ); % (ode)met
% viral ODE simulation
odeOpt = odeset('AbsTol',1e-3);
odeSolver = organism.methods.viralODEsolverHandle;
[~, cViral] = odeSolver( viralODEfn, tViralLim, cViralInitial,odeOpt, organism.modelInfected.viralVarargin, viralODElimits ); % sec, (ode)met
% convert to flux over time

fViral = ((cViral(end,:)' - cViralInitial) / tStep) * organism.modelInfected.ode2fbaMetUnits ; % ((ode)met / hour)  to ((fba)met / hour)


% constrain host fluxes
organism.simulation.fluxBoundsLower(organism.model.viralRxnODETranslation(:,1),iSim) = fViral(organism.model.viralRxnODETranslation(:,2));
organism.simulation.fluxBoundsUpper(organism.model.viralRxnODETranslation(:,1),iSim) = fViral(organism.model.viralRxnODETranslation(:,2));
% store new virus concentrations
% store in output for each FBA step
organism.simulation.viralConcentrations( : , iSim ) = cViral(end,:)' ; %
% also store in output with viral substeps
if size( cViral(2:end,:)', 2) == 0
    % replace with zero if ode solver fails.
    cViral(2 : organism.simulation.nViralSubSteps + 1 , 1:organism.model.nViralODEs) = zeros(  organism.simulation.nViralSubSteps , organism.model.nViralODEs );
end
organism.simulation.viralSubSteps( :, organism.simulation.nViralSubSteps * (iSim - 2) + 2 : organism.simulation.nViralSubSteps * (iSim-1) + 1 ) = cViral(2:end,:)';
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

organismSimulation = organism.simulation;