optEnvelope

src/design/optEnvelope/addEnv.m

@@ -0,0 +1,111 @@
+function line = addEnv(origModel, biomass, desiredProduct, varargin)
+% addEnv adds envelope to figure
+% Algorithm is able to knock out genes as well as reactions to produce
+% production envelope
+%
+% USAGE:
+%   line = addEnv(origModel, biomass, desiredProduct, 'KnockOuts', knockouts, 'colour', colour, 'prodMol', prodMol, 'subUptake', subUptake, 'molarSum', molarSum)
+%
+% INPUTS:
+%   origModel         COBRA model structure [struct]
+%   biomass           Reaction name of biomass [char]
+%   desiredProduct    Reaction name of desired product [char]
+%
+% OPTIONAL INPUTS:
+%   KnockOuts         (opt) List of knockouts for production envelope [cell array] (default: {})
+%   colour            (opt) Short name for colour of line to plot [anything that can be colour in matlab] (default: 'r' (red))
+%   prodMol           (opt) Molar mass of target product for yield plot [double]
+%   subUptake         (opt) Uptake of substrate for yield plot [double]
+%   molarSum          (opt) Molar mass of substrate for yield plot [double]
+%
+% OUTPUTS
+%   line              Line data for plot function
+%
+% NOTES
+%   Sometimes last point of envelope drops to zero (might be rounding error)
+%   but this function connects last points of lines so the graph creates
+%   continuous line.
+%   This algorithm only adds graph. It does not change labels.
+%
+% EXAMPLE:
+%   line = addEnv(model, 'BIOMASS_Ecoli', 'EX_ac_e', {'GHMT2r','GND','SUCCt2r','SUCD4','AKGt2r','GLUt2r'}, 'm')
+%
+% AUTHORS:
+%   Created by Kristaps Berzins    31/10/2022
+%   Modified by Kristaps Berzins   30/09/2024
+
+parser = inputParser();
+parser.addRequired('model', @(x) isstruct(x) && isfield(x, 'S') && isfield(origModel, 'rxns')...
+    && isfield(origModel, 'mets') && isfield(origModel, 'lb') && isfield(origModel, 'ub') && isfield(origModel, 'b')...
+    && isfield(origModel, 'c'))
+parser.addRequired('biomass', @(x) any(validatestring(x, origModel.rxns)))
+parser.addRequired('desiredProduct', @(x) any(validatestring(x, origModel.rxns)))
+parser.addOptional('KnockOuts', {}, @(x) iscell(x) && ismatrix(x))
+parser.addOptional('colour', 'r', @(x) any(validatecolor(x)))
+parser.addOptional('prodMol', [], @(x) isnumeric(x))
+parser.addOptional('subUptake', 10, @(x) isnumeric(x))
+parser.addOptional('molarSum', 180, @(x) isnumeric(x))
+
+parser.parse(origModel, biomass, desiredProduct, varargin{:});
+origModel = parser.Results.model;
+biomass = parser.Results.biomass;
+desiredProduct = parser.Results.desiredProduct;
+KnockOuts = parser.Results.KnockOuts;
+colour = parser.Results.colour;
+prodMol = parser.Results.prodMol;
+subUptake = parser.Results.subUptake;
+molarSum = parser.Results.molarSum;
+
+if isempty(prodMol)
+    prodMolIs = false;
+else
+    prodMolIs = true;
+end
+
+model = origModel;
+
+if any(ismember(model.rxns, KnockOuts))
+    rxns = ismember(model.rxns, KnockOuts);
+    model.ub(rxns) = 0;
+    model.lb(rxns) = 0;
+elseif any(ismember(model.genes, KnockOuts))
+    model = buildRxnGeneMat(model);
+    [model, ~, ~] = deleteModelGenes(model, KnockOuts);
+%elseif %Enzymes
+end
+
+solMin = optimizeCbModel(model, 'min');
+solMax = optimizeCbModel(model, 'max');
+controlFlux1 = linspace(solMin.f, solMax.f, 100)';
+if nnz(controlFlux1) == 0
+    return;
+end
+model = changeObjective(model, desiredProduct);
+
+for i = 1:numel(controlFlux1)
+    model = changeRxnBounds(model, biomass, controlFlux1(i), 'b');
+    s = optimizeCbModel(model, 'min'); Min1(i, 1) = s.f;
+    if s.stat == 0
+        model = changeRxnBounds(model, biomass, controlFlux1(i) - 0.0001 * controlFlux1(i), 'b');
+        s = optimizeCbModel(model, 'min'); Min1(i, 1) = s.f;
+        s = optimizeCbModel(model, 'max'); Max1(i, 1) = s.f;
+    end
+    s = optimizeCbModel(model, 'max'); Max1(i, 1) = s.f;
+    if s.stat == 0
+        model = changeRxnBounds(model, biomass, controlFlux1(i) - 0.0001 * controlFlux1(i), 'b');
+        s= optimizeCbModel(model,'min');Min1(i,1)=s.f;
+        s= optimizeCbModel(model,'max');Max1(i,1)=s.f;
+    end
+end
+
+if prodMolIs
+    controlFlux1 = controlFlux1 / subUptake * 1000 / molarSum;
+    Max1 = Max1 / molarSum * prodMol / subUptake;
+    Min1 = Min1 / molarSum * prodMol / subUptake;
+end
+
+hold on
+line = plot(controlFlux1, Max1, 'color', colour, 'LineWidth', 2);
+plot(controlFlux1, Min1, 'color', colour, 'LineWidth', 2)
+hold off
+

src/design/optEnvelope/milpOEReinserts.m

@@ -0,0 +1,188 @@
+function [knockouts] = milpOEReinserts(model, data, K, minP, numKO, toDel, timeLimit, printLevel)
+% This function is creating and calculating MILP to get certain amount
+% of knockouts to achieve best production envelope
+%
+% USAGE:
+%   [knockouts] = milpOEReinserts(model, data, K, minP, numKO, toDel, timeLimit, printLevel)
+%
+% INPUTS:
+%   model       COBRA model structure in irreversible form [struct]
+%   data        Struct with information about:
+%                   * mainActive    List of active reactions for main envelope [cell array]
+%   K           List of reactions that cannot be selected for knockout (reaction IDs) [double array]
+%   minP        Struct with information about biomass and desired product.
+%               	* bioID         ID of biomass [double]
+%                   * proID         ID of desired product[double]
+%   numKO       Number of knockouts to achieve [double]
+%
+% OPTIONAL INPUTS:
+%   toDel       Numeric variable that shows what to delete:
+%                   0: reactions
+%                   1: genes
+%                   2: enzymes
+%   timeLimit   Time limit for gurobi optimization (in seconds) [double]
+%   printLevel  Print level for gurobi optimization [double]
+%
+% OUTPUTS:
+%   knockouts   List of reactions that when removed gives optimal envelope
+%
+% AUTHORS:
+%   created by Kristaps Berzins 31/10/2022
+%
+% NOTES:
+%   This function is not designed for stand-alone use. Should be used by
+%   using optEnvelope.m with set numKO parameter
+
+if nargin < 6
+    toDel = 0;
+end
+if nargin < 7
+    timeLimit = inf;
+end
+if nargin < 8
+    printLevel = 0;
+end
+%%
+switch toDel
+    case 0
+        rxns = findRxnIDs(model, data.mainActive);
+        rxns = ismember(model.rxns,[model.rxns(K);model.rxns(rxns)]);
+        reactions = model.rxns(~rxns);
+
+        % MILP
+
+        [r,c] = size(model.S);
+
+        model.C_chemical=zeros(c,1);
+        model.C_chemical(minP.proID)=1;
+
+        yInd=[];
+        notYInd=[];
+        reactionIDs = findRxnIDs(model, reactions);
+
+        for i=1:c
+            if find(reactionIDs == i)
+                yInd(end+1)=i;
+                continue;
+            else
+                notYInd(end+1)=i;
+            end
+        end
+
+        notYInd=sort(notYInd)';
+        yInd=sort(yInd)';
+
+        ySize=size(yInd,1);
+        I = eye(c);
+        A=[model.S;-model.S;I(notYInd,:);-I(notYInd,:); I(yInd,:); -I(yInd,:)];
+
+        [aSizeRow, vSize]=size(A);
+        reinserts=zeros(c,1);
+        reinserts(yInd)=1;
+
+        Ay1=diag(reinserts); %ub
+        Ay1=Ay1*diag(model.ub);
+
+        Ay2=diag(reinserts); %lb
+        Ay2=Ay2*diag(model.lb);
+
+        z1=find(Ay1);
+        z2=find(Ay2);
+        zSize=size([z1;z2],1);
+
+        Ay= [zeros(2*r+2*(vSize-ySize),ySize);
+            -Ay1(yInd,yInd);
+            Ay2(yInd,yInd);
+            ];
+
+        B = [zeros(2*r,1);
+            model.ub(notYInd);
+            -model.lb(notYInd);
+            zeros(2*(ySize),1);
+            ];
+
+        C=model.c;
+
+        [A_w,Ay_w ,B_w,C_w, ~, ~, wSize] = seperateTransposeJoinOE(-A,-Ay,-B,-C,ySize,1,c,1000,zSize); %max(model.ub) 1000
+
+        awSizeRow=size(A_w,1);
+
+        Cjoined=[model.C_chemical; zeros(wSize+zSize,1); zeros(ySize,1)];
+
+        A2=-[
+            -C', -C_w';
+            A, sparse(aSizeRow, wSize+zSize);
+            sparse(awSizeRow, vSize), A_w];
+
+        Ay2=-[
+            zeros(1, ySize);
+            Ay;
+            Ay_w];
+
+        C2=Cjoined(1:vSize+wSize+zSize,:);
+
+        B2=-[
+            0;
+            B;
+            B_w;
+            ];
+
+        z3=find(Ay2);
+        zSizeOptKnock2=size(z3,1);
+
+        [A2_w,Ay2_w,B2_w,C2_w,lb2_w,ub2_w,uSize]=seperateTransposeJoinOE(A2,Ay2,B2,C2,ySize,1,vSize+wSize+zSize,1000,zSizeOptKnock2); %max(model.ub) 1000
+
+        [A2_wRow, A2_wCol]=size(A2_w);
+        [ARow, ACol]=size(A);
+
+        A3=[
+            A2_w, sparse(A2_wRow,ACol), Ay2_w;                      %dual constraints
+            zeros(1,uSize+zSizeOptKnock2+vSize),  -ones(1,ySize);   %y sum constraints
+            sparse(ARow,uSize+zSizeOptKnock2), A, Ay;               %feasibility conatraint
+            ];
+
+        B3=[
+            B2_w;
+            numKO-ySize;
+            B;
+            ];
+
+        C3=[C2_w;
+            zeros(ACol,1);
+            zeros(ySize,1);
+            ];
+
+        tmpL=lb2_w(1:uSize+zSizeOptKnock2);
+        tmpH=ub2_w(1:uSize+zSizeOptKnock2);
+        ysUpperBound=ones(ySize,1);
+        lb3=[tmpL; model.lb; zeros(ySize,1)];
+        ub3=[tmpH; model.ub; ysUpperBound];
+        intVars=A2_wCol+ACol+1:A2_wCol+ACol+ySize;
+
+        MILP.c = C3;
+        MILP.osense = -1; % max
+        MILP.A = sparse(A3);
+        MILP.b = B3;
+        MILP.lb = lb3; MILP.ub = ub3;
+        MILP.x0 = [];
+        MILP.vartype = char(ones(1,length(C3)).*double('C'));
+        MILP.vartype(intVars) = 'B';
+        MILP.csense = char(ones(1,length(B3)).*double('L'));
+
+        milpSol = solveCobraMILP(MILP, 'timeLimit', timeLimit, 'printLevel', printLevel);
+
+        %% add only found reactions
+
+        rxns = reactions(milpSol.int == 1);
+        knockouts = reactions(milpSol.int ~= 1);
+        numel(rxns)
+        numel(knockouts)
+        %debug
+        [irrev,~,~,~] = convertToIrreversible(model);
+        addEnv(irrev, irrev.rxns(minP.bioID), irrev.rxns(minP.proID), knockouts, 'r');
+        print('debug');
+    case 1
+        %Genes
+    case 2
+        %Enzymes
+end