Skip to content

Commit

Permalink
Merge pull request #2330 from Kretaks/optEnvelope
Browse files Browse the repository at this point in the history 
optEnvelope
  • Loading branch information
@rmtfleming
rmtfleming authored Nov 7, 2024
2 parents c8acdba + 6cdc246 commit fa2adac
Show file tree
Hide file tree
Showing 8 changed files with 1,042 additions and 0 deletions.
111 changes: 111 additions & 0 deletions src/design/optEnvelope/addEnv.m
View file
Original file line number Diff line number Diff line change
@@ -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

188 changes: 188 additions & 0 deletions src/design/optEnvelope/milpOEReinserts.m
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit fa2adac

Please sign in to comment.