optimizeRegModel

PURPOSE ^

optimizeRegModel - finds the steady state solution of a model with

SYNOPSIS ^

function [FBAsols,DRgenes,constrainedRxns,cycleStart,states] = optimizeRegModel(model,initialRegState)

DESCRIPTION ^

 optimizeRegModel - finds the steady state solution of a model with
 Boolean regulatory constraints

 [FBAsols,DRgenes,constrainedRxns,cycleStart,states] = optimizeRegModel(model,initialRegState)

 model             a regulatory COBRA model
 initialRegState   the initial state of the regulatory network as a 
                   Boolean vector (opt, default = all false)

 FBAsols           all of the FBA solutions at the steady state (or stable
                   cycle) of the regulatory network 
 DRgenes           the genes that are OFF for every FBA solution
 constrainedRxns   the reactions that are OFF for every FBA solution
 cycleStart        the number of iterations before the regulatory network
                   reaches the steady state or cycle
 states            the state of the regulatory network at every iteration
                   calculated

 Jeff Orth 8/19/08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [FBAsols,DRgenes,constrainedRxns,cycleStart,states] = optimizeRegModel(model,initialRegState)
0002 % optimizeRegModel - finds the steady state solution of a model with
0003 % Boolean regulatory constraints
0004 %
0005 % [FBAsols,DRgenes,constrainedRxns,cycleStart,states] = optimizeRegModel(model,initialRegState)
0006 %
0007 % model             a regulatory COBRA model
0008 % initialRegState   the initial state of the regulatory network as a
0009 %                   Boolean vector (opt, default = all false)
0010 %
0011 % FBAsols           all of the FBA solutions at the steady state (or stable
0012 %                   cycle) of the regulatory network
0013 % DRgenes           the genes that are OFF for every FBA solution
0014 % constrainedRxns   the reactions that are OFF for every FBA solution
0015 % cycleStart        the number of iterations before the regulatory network
0016 %                   reaches the steady state or cycle
0017 % states            the state of the regulatory network at every iteration
0018 %                   calculated
0019 %
0020 % Jeff Orth 8/19/08
0021 
0022 if nargin < 2
0023     rFBAsol1 = false.*ones(length(model.regulatoryGenes),1); % initial state
0024     inputs1state = false.*ones(length(model.regulatoryInputs1),1);
0025     inputs2state = false.*ones(length(model.regulatoryInputs2),1);
0026     state1 = false.*ones(length(model.regulatoryGenes)+length(model.regulatoryInputs1)+length(model.regulatoryInputs2),1); %vector for entire state
0027 else
0028     if length(initialRegState) ~= (length(model.regulatoryGenes)+length(model.regulatoryInputs1)+length(model.regulatoryInputs2))
0029         error('initialRegState is invalid length');
0030     end
0031     state1 = initialRegState;
0032     rFBAsol1 = state1(1:length(model.regulatoryGenes));
0033     inputs1state = state1((length(model.regulatoryGenes)+1):(length(model.regulatoryGenes)+length(model.regulatoryInputs1)));
0034     inputs2state = state1((length(model.regulatoryGenes)+length(model.regulatoryInputs1)+1):(length(model.regulatoryGenes)+length(model.regulatoryInputs1)+length(model.regulatoryInputs2)));
0035 end
0036     
0037 rFBAsols = {rFBAsol1}; % matrices of all solutions
0038 inputs1states = {inputs1state};
0039 inputs2states = {inputs2state};
0040 states = {state1};
0041 
0042 [rFBAsol2,finalInputs1States,finalInputs2States] = solveBooleanRegModel(model,rFBAsol1,inputs1state,inputs2state); %get first solution
0043 rFBAsols{2} = rFBAsol2; % add solution to array of solutions
0044 inputs1states{2} = finalInputs1States;
0045 inputs2states{2} = finalInputs2States;
0046 states{2} = [rFBAsol2;finalInputs1States;finalInputs2States];
0047 
0048 % continue solving until steady state
0049 cycleReached = false;
0050 while ~cycleReached % while current solution ~= any previous solutions
0051     rFBAsol1 = rFBAsols{length(rFBAsols)};
0052     inputs1state = inputs1states{length(inputs1states)};
0053     inputs2state = inputs2states{length(inputs2states)};
0054     [rFBAsol2,finalInputs1States,finalInputs2States] = solveBooleanRegModel(model,rFBAsol1,inputs1state,inputs2state);
0055     rFBAsols{length(rFBAsols)+1} = rFBAsol2; % add solution to array of solutions
0056     inputs1states{length(inputs1states)+1} = finalInputs1States;
0057     inputs2states{length(inputs2states)+1} = finalInputs2States;
0058     states{length(states)+1} = [rFBAsol2;finalInputs1States;finalInputs2States];
0059     %check if this state has been reached before
0060     cycleStart = [];
0061     for i = 1:(length(states)-1)
0062         if all(states{length(states)} == states{i})
0063             cycleStart = i;
0064         end
0065     end
0066     if any(cycleStart)
0067         cycleReached = true;
0068     end
0069 end
0070 
0071 % remove downregulated genes and compute growth by FBA
0072 FBAsols = {};
0073 DRgenes = {};
0074 constrainedRxns = {};
0075 k = 0;
0076 for i = cycleStart:(length(states)-1)
0077     k = k+1;
0078     genes = {};
0079         for j = 1:length(model.regulatoryGenes)
0080             if rFBAsols{i}(j) == false 
0081                 genes{length(genes)+1,1} = model.regulatoryGenes{j};
0082             end
0083         end
0084     genes = intersect(model.genes,genes); % remove genes not associated with rxns
0085     [modelDR,he,rxns] = deleteModelGenes(model,genes); % set rxns to 0
0086     fbasol = optimizeCbModel(modelDR,'max',true);
0087     FBAsols{k} = fbasol;
0088     DRgenes{k} = genes;
0089     constrainedRxns{k} = rxns;
0090 end
0091 
0092 
0093 
0094

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