% 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: set the current regulatory state by evaluation of rules for genes and TFs
% 
%
% Inputs: 
%
% Outputs: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function [tfActivityCurrent,geneActivityCurrent,geneProductActivityCurrent] = ...
    FBA_EB_setCurrentRegulation(tStep,iSim,organism,organismMedia,simulationStimuli)
                            
% this function determines the current activity of regulatory elements: tfs, genes, and
% gene products

simulationGenes = organism.simulation.genes;
simulationGeneProducts = organism.simulation.geneProducts;
simulationTFs = organism.simulation.tfs;
simulationFluxes = organism.simulation.fluxes;
currentMedia = organismMedia; % this is a reference to the previous timestep, iSim - 1, resources allocated to this organism
currentStimuli = simulationStimuli; % this is a reference to the previous timestep, iSim - 1
organismModelStruct = organism.model;

% initialize output variables to fill
tfActivityCurrent = false(organismModelStruct.nTFs,1);
geneActivityCurrent = false(organismModelStruct.nGene,1);
geneProductActivityCurrent = false(organismModelStruct.nGene,1);

% first extract relevent data for rule evaluation from simulation states
% Note, currently all these rules reference the state of the previous simulation, to avoid
% self reference durring single timestep.
% Note also, the names of the following variables are required to match the rule strings
% used later in eval( ), if change is desired address in rule parser also.
geneProductActivity = simulationGeneProducts(:,iSim - 1);
geneActivity = simulationGenes(:,iSim-1);
currentFlux = simulationFluxes(:,iSim - 1);
tfActivity = simulationTFs(:,iSim - 1);

% for the purpose of regulation, round very small media/fluxes to zero
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
zeroThresh = 0.0001;

currentFlux = currentFlux .* ( abs( currentFlux ) > zeroThresh );
currentMedia = currentMedia .* ( abs( currentMedia ) > zeroThresh );

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

% Evaluate gene activity rules
for iGeneAct = 1 : organismModelStruct.nGene
    geneActivityCurrent(iGeneAct) = eval(organismModelStruct.geneRegulationRulesEvalStrings{iGeneAct});
end % for iGeneAct

% Evaluate TF activity rules
for iTFAct = 1 : organismModelStruct.nTFs
    tfActivityCurrent(iTFAct) = eval(organismModelStruct.tfRegulationRulesEvalStrings{iTFAct});
end % for iTFAct

% Apply delay to past gene activity to determine current gene product activity
for iGeneProdAct = 1 : organismModelStruct.nGene
    delayTimesteps = floor(organismModelStruct.activityDelayList(iGeneProdAct)/tStep);
    % check if delay referrences a time before simulation start
    if delayTimesteps < iSim
        % delay references time durring simulation
        geneProductActivityCurrent(iGeneProdAct) = simulationGenes(iGeneProdAct, (iSim - delayTimesteps) );
    else
        % delay references time before simulation initialized, default to true
        geneProductActivityCurrent(iGeneProdAct) = simulationGenes(iGeneProdAct, 1 );
    end % if delayTimesteps < iSim
    
end % for iGeneProdAct
end % FBA_EB_setCurrentRegulation
