solveBooleanRegModel

PURPOSE ^

solveBooleanRegModel - determines the next state of the regulatory

SYNOPSIS ^

function [finalState,finalInputs1States,finalInputs2States] = solveBooleanRegModel(model,initialState,inputs1States,inputs2States)

DESCRIPTION ^

 solveBooleanRegModel - determines the next state of the regulatory
 network based on the current state. Called by optimizeRegModel and
 dynamicRFBA

 [finalState,finalInputs1States,finalInputs2States] =
 solveBooleanRegModel(model,initialState,inputs1States,inputs2States)

 model                 a regulatory COBRA model
 initialState          initial state of regulatory network
 inputs1States         initial state of type 1 inputs (metabolites)
 inputs2States         initial state of type 2 inputs (reactions)

 finalState            final state of regulatory network
 finalInputs1States    final state of type 1 inputs
 finalInputs2States    final state of type 2 inputs

 Jeff Orth  7/24/08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [finalState,finalInputs1States,finalInputs2States] = solveBooleanRegModel(model,initialState,inputs1States,inputs2States)
0002 % solveBooleanRegModel - determines the next state of the regulatory
0003 % network based on the current state. Called by optimizeRegModel and
0004 % dynamicRFBA
0005 %
0006 % [finalState,finalInputs1States,finalInputs2States] =
0007 % solveBooleanRegModel(model,initialState,inputs1States,inputs2States)
0008 %
0009 % model                 a regulatory COBRA model
0010 % initialState          initial state of regulatory network
0011 % inputs1States         initial state of type 1 inputs (metabolites)
0012 % inputs2States         initial state of type 2 inputs (reactions)
0013 %
0014 % finalState            final state of regulatory network
0015 % finalInputs1States    final state of type 1 inputs
0016 % finalInputs2States    final state of type 2 inputs
0017 %
0018 % Jeff Orth  7/24/08
0019 
0020 
0021 % determine state of inputs
0022 
0023 % determine external metabolite levels from exchange rxn bounds (maybe change this later)
0024 finalInputs1States = [];
0025 [selExc,selUpt] = findExcRxns(model); %get all exchange rxns
0026 for i = 1:length(model.regulatoryInputs1)
0027     met = model.regulatoryInputs1{i};
0028     fullS = full(model.S);
0029     rxnID = intersect(find(fullS(findMetIDs(model,met),:)),find(selExc));
0030     if model.lb(rxnID) < 0
0031         finalInputs1States(i,1) = true;
0032     else
0033         finalInputs1States(i,1) = false;
0034     end
0035 end    
0036 
0037 % apply initialState to model, get rxn fluxes
0038 drGenes = {};
0039 for i = 1:length(model.regulatoryGenes)
0040     if initialState(i) == false 
0041         drGenes{length(drGenes)+1} = model.regulatoryGenes{i};
0042     end
0043 end
0044 drGenes = intersect(model.genes,drGenes); % remove genes not associated with rxns
0045 modelDR = deleteModelGenes(model,drGenes); % set rxns to 0
0046 fbasolDR = optimizeCbModel(modelDR,'max',true);
0047 
0048 % if rxn flux = 0, set state to false
0049 finalInputs2States = [];
0050 if ~any(fbasolDR.x)
0051     finalInputs2States = false.*ones(length(model.regulatoryInputs2),1);
0052 else
0053     for i = 1:length(model.regulatoryInputs2)
0054         rxnFlux = fbasolDR.x(findRxnIDs(model,model.regulatoryInputs2{i}));
0055         if rxnFlux == 0
0056             finalInputs2States(i,1) = false;
0057         else
0058             finalInputs2States(i,1) = true;
0059         end
0060     end
0061 end
0062 
0063 % determine state of genes
0064 
0065 ruleList = parseRegulatoryRules(model); %get the set of rules in a form that can be evaluated
0066 
0067 geneStates = [];
0068 for i = 1:length(model.regulatoryRules)
0069     geneStates(i,1) = eval(ruleList{i});
0070 end
0071 
0072 finalState = geneStates;
0073     
0074 
0075 %% parseRegulatoryRules
0076 function ruleList = parseRegulatoryRules(model)
0077 
0078 ruleList = cell(length(model.regulatoryRules),1); %preallocate array
0079 
0080 for i = 1:length(model.regulatoryRules)
0081     fields = splitString(model.regulatoryRules{i},'[\s~|&()]');
0082     newFields = fields;
0083     for j = 1:length(fields) %iterate through words and replace
0084         word = fields{j};
0085         if strcmp(word,'true') 
0086             newFields{j} = 'true';
0087         elseif strcmp(word,'false')
0088             newFields{j} = 'false';
0089         else
0090             if any(strcmp(word,model.regulatoryGenes))
0091                 index = find(strcmp(word,model.regulatoryGenes));
0092                 newFields{j} = ['initialState(',num2str(index),')'];
0093             elseif any(strcmp(word,model.regulatoryInputs1))
0094                 index = find(strcmp(word,model.regulatoryInputs1));
0095                 newFields{j} = ['inputs1States(',num2str(index),')'];
0096             elseif any(strcmp(word,model.regulatoryInputs2))
0097                 index = find(strcmp(word,model.regulatoryInputs2));
0098                 newFields{j} = ['inputs2States(',num2str(index),')'];
0099             else
0100                 newFields{j} = ''; %no match, delete the invalid word
0101             end
0102         end
0103     end
0104     newRule = model.regulatoryRules{i};
0105     for j = 1:length(fields)
0106         newRule = strrep(newRule,fields{j},newFields{j});
0107     end
0108     ruleList{i} = newRule;
0109 end
0110     
0111

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003