diff --git a/tutorials/optKnock/tutorial_optKnock.html b/tutorials/optKnock/tutorial_optKnock.html new file mode 100644 index 0000000000..377c61bdb3 --- /dev/null +++ b/tutorials/optKnock/tutorial_optKnock.html @@ -0,0 +1,557 @@ + +OptKnock Tutorial

OptKnock Tutorial

Author: Sebastián N. Mendoza, Center for Mathematical Modeling, University of Chile. snmendoz@uc.cl

Reviewer(s): Thomas Pfau, Sylvian Arreckx, Laurent Heirendt & Ronan Fleming

INTRODUCTION:

In this tutorial we will run optKnock. For a detailed description of the procedure, please see [1]. Briefly, the problem is to find a set of reactions of size "numDel" such that when these reactions are deleted from the network, the mutant created will produce a particular target of interest in a higher rate than the wild-type strain.

For example, imagine that we would like to increase the production of succinate or lactate in Escherichia coli. Which are the knock-outs needed to increase the production of these products? We will approach these problems in this tutorial.

MATERIALS

EQUIPMENT

  1. MATLAB
  2. A solver for Mixed Integer Linear Programming (MILP) problems. For example, Gurobi.

EQUIPMENT SETUP

Use changeCobraSolver to choose the solver for MILP problems.

PROCEDURE

The proceduce consists on the following steps

1) Define contraints

2) Define a set of reactions to search knockouts. Only reactions in this set will be deleted.

3) Define the number of reactions to be deleter, the target reaction and some constraints to be accomplish

4) Run optKnock. TIMING: This task should take from a few seconds to a few hours depending on the size of your reconstruction

We verify that cobratoolbox has been initialized and that the solver has been set.

global TUTORIAL_INIT_CB;

if ~isempty(TUTORIAL_INIT_CB) && TUTORIAL_INIT_CB==1

initCobraToolbox

changeCobraSolver('gurobi','all');

end

fullPath = which('optKnockTutorial');

folder = fileparts(fullPath);

currectDirectory = pwd;

cd(folder);

We load the model of E. coli [2]

model = readCbModel('iJO1366.mat')

biomass = 'BIOMASS_Ec_iJO1366_core_53p95M';

We define the maximum number of solutions to find

threshold = 5;

First, we define the set for reactions which could be deleted from the network. Reactions not in this list are not going to be deleted.

selectedRxnList = {'GLCabcpp'; 'GLCptspp'; 'HEX1'; 'PGI'; 'PFK'; 'FBA'; 'TPI'; 'GAPD'; 'PGK'; 'PGM'; 'ENO'; 'PYK'; 'LDH_D'; 'PFL'; 'ALCD2x'; 'PTAr'; 'ACKr'; 'G6PDH2r'; 'PGL'; 'GND'; 'RPI'; 'RPE'; 'TKT1'; 'TALA'; 'TKT2'; 'FUM'; 'FRD2'; 'SUCOAS'; 'AKGDH'; 'ACONTa'; 'ACONTb'; 'ICDHyr'; 'CS'; 'MDH'; 'MDH2'; 'MDH3'; 'ACALD'};

Then, we define some constraints

% prespecified amount of glucose uptake 10 mmol/grDW*hr

model = changeRxnBounds(model, 'EX_glc__D_e', -10, 'b');

% Unconstrained uptake routes for inorganic phosphate, sulfate and

% ammonia

model = changeRxnBounds(model, 'EX_o2_e', 0, 'l');

model = changeRxnBounds(model, 'EX_pi_e', -1000, 'l');

model = changeRxnBounds(model, 'EX_so4_e', -1000, 'l');

model = changeRxnBounds(model, 'EX_nh4_e', -1000, 'l');

% The optimization step could opt for or against the phosphotransferase

% system, glucokinase, or both mechanisms for the uptake of glucose

model = changeRxnBounds(model, 'GLCabcpp', -1000, 'l');

model = changeRxnBounds(model, 'GLCptspp', -1000, 'l');

model = changeRxnBounds(model, 'GLCabcpp', 1000, 'u');

model = changeRxnBounds(model, 'GLCptspp', 1000, 'u');

model = changeRxnBounds(model, 'GLCt2pp', 0, 'b');

% Secretion routes for acetate, carbon dioxide, ethanol, formate, lactate

% and succinate are enabled

model = changeRxnBounds(model, 'EX_ac_e', 1000, 'u');

model = changeRxnBounds(model, 'EX_co2_e', 1000, 'u');

model = changeRxnBounds(model, 'EX_etoh_e', 1000, 'u');

model = changeRxnBounds(model, 'EX_for_e', 1000, 'u');

model = changeRxnBounds(model, 'EX_lac__D_e', 1000, 'u');

model = changeRxnBounds(model, 'EX_succ_e', 1000, 'u');

Then, we calculates the production of metabolites before running optKnock

% determine succinate production and growth rate before optimizacion

fbaWT = optimizeCbModel(model);

succFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_succ_e'));

etohFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_etoh_e'));

formFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_for_e'));

lactFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_lac__D_e'));

acetFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_ac_e'));

growthRateWT = fbaWT.f;

fprintf('The production of succinate before optimization is %.1f \n', succFluxWT);

The production of succinate before optimization is 0.1

fprintf('The growth rate before optimization is %.1f \n', growthRateWT);

The growth rate before optimization is 0.2

fprintf('The production of other products such as ethanol, formate, lactate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxWT, formFluxWT, lactFluxWT, acetFluxWT);

The production of other products such as ethanol, formate, lactate and acetate are 8.1, 17.3, 0.0 and 8.2, respectively.

I) SUCCINATE OVERPRODUCTION

EXAMPLE 1: finding optKnock reactions sets of large 2 for increasing production of succinate

fprintf('\n...EXAMPLE 1: Finding optKnock sets of large 2 or less...\n\n')

...EXAMPLE 1: Finding optKnock sets of large 2 or less...

% Set optKnock options

% The exchange of succinate will be the objective of the outer problem

options = struct('targetRxn', 'EX_succ_e', 'numDel', 2);

% We will impose that biomass be at least 50% of the biomass of wild-type

constrOpt = struct('rxnList', {{biomass}},'values', 0.5*fbaWT.f, 'sense', 'G');

% We will try to find 10 optKnock sets of a maximun length of 2

previousSolutions = cell(10, 1);

contPreviousSolutions = 1;

nIter = 1;

while nIter < threshold

fprintf('...Performing optKnock analysis...\n')

if isempty(previousSolutions{1})

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt);

else

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions, 1);

end

% determine succinate production and growth rate after optimizacion

succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e'));

growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns, biomass));

etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e'));

formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e'));

lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e'));

acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e'));

setM1 = optKnockSol.rxnList;

if ~isempty(setM1)

previousSolutions{contPreviousSolutions} = setM1;

contPreviousSolutions = contPreviousSolutions + 1;

%printing results

fprintf('optKnock found a optKnock set of large %d composed by ', length(setM1));

for j = 1:length(setM1)

if j == 1

fprintf('%s', setM1{j});

elseif j == length(setM1)

fprintf(' and %s', setM1{j});

else

fprintf(', %s', setM1{j});

end

end

fprintf('\n');

fprintf('The production of succinate after optimization is %.2f \n', succFluxM1);

fprintf('The growth rate after optimization is %.2f \n', growthRateM1);

fprintf('The production of other products such as ethanol, formate, lactate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, lactFluxM1, acetFluxM1);

fprintf('...Performing coupling analysis...\n');

[type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_succ_e');

fprintf('The solution is of type: %s\n', type);

fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth);

fprintf('The maximun and minimun production of succinate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd);

if strcmp(type, 'growth coupled')

singleProductionEnvelope(model, setM1, 'EX_succ_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['succ_ex1_' num2str(nIter)], 'outputFolder', 'OptKnockResults');

end,

else

if nIter == 1

fprintf('optKnock was not able to found an optKnock set\n');

else

fprintf('optKnock was not able to found additional optKnock sets\n');

end

break;

end

nIter = nIter + 1;

end

...Performing optKnock analysis...
optKnock found a optKnock set of large 2 composed by
PFL
and TKT2
The production of succinate after optimization is 1.68
The growth rate after optimization is 0.22
The production of other products such as ethanol, formate, lactate and acetate are 14.9, 0.0, -0.0 and -0.0, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.18
The maximun and minimun production of succinate given the optKnock set is 0.06 and 1.41, respectively
...Performing optKnock analysis...
MILP problem with 8398 constraints 63 integer variables and 8370 continuous variables
Optimize a model with 8398 rows, 8370 columns and 34671 nonzeros
Variable types: 8307 continuous, 63 integer (63 binary)
Coefficient statistics:
  Matrix range     [2e-06, 1e+03]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+03]
  RHS range        [1e-02, 1e+03]
Presolve removed 5079 rows and 3321 columns
Presolve time: 0.28s
Presolved: 3319 rows, 5049 columns, 21617 nonzeros
Variable types: 5012 continuous, 37 integer (37 binary)

Root relaxation: objective 1.230321e+01, 3505 iterations, 0.79 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   12.30321    0    2          -   12.30321      -     -    2s
H    0     0                       0.0915893   12.30321      -     -    2s
     0     0   12.30321    0    2    0.09159   12.30321      -     -    2s
     0     2   12.30321    0    2    0.09159   12.30321      -     -    2s
    81    71    5.47336   24    4    0.09159   12.30321      -  64.2    5s
*   91    62              22       0.1718801   12.30321  7058%  73.0    5s
*  167    71              25       0.2411614   12.30321  5002%  81.0    8s
H  170    71                       0.2864668   12.30321  4195%  79.6    8s
   194    77   12.30321   19    6    0.28647   12.30321  4195%  82.5   10s
*  242    89              23       1.2019176   12.30321   924%  80.7   10s
*  292    91              23       1.4425937   12.29845   753%  81.0   12s
   339    37 infeasible   17         1.44259   10.09393   600%  84.7   15s
*  361    22              21       2.7454999    8.37613   205%  84.5   15s

Explored 398 nodes (39947 simplex iterations) in 17.01 seconds
Thread count was 4 (of 4 available processors)

Solution count 7: 2.7455 1.44259 1.20192 ... 0.0915893
Pool objective bound 2.7455

Optimal solution found (tolerance 1.00e-12)
Warning: max constraint violation (4.4603e-09) exceeds tolerance
Best objective 2.745499915504e+00, best bound 2.745499915504e+00, gap 0.0000%
optKnock found a optKnock set of large 2 composed by
PFL
and RPI
The production of succinate after optimization is 2.75
The growth rate after optimization is 0.14
The production of other products such as ethanol, formate, lactate and acetate are 0.0, 0.0, 9.4 and -0.0, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.12
The maximun and minimun production of succinate given the optKnock set is 0.04 and 2.31, respectively
...Performing optKnock analysis...
MILP problem with 8399 constraints 63 integer variables and 8370 continuous variables
Optimize a model with 8399 rows, 8370 columns and 34674 nonzeros
Variable types: 8307 continuous, 63 integer (63 binary)
Coefficient statistics:
  Matrix range     [2e-06, 1e+03]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+03]
  RHS range        [1e-02, 1e+03]
Presolve removed 5079 rows and 3321 columns
Presolve time: 0.27s
Presolved: 3320 rows, 5049 columns, 21619 nonzeros
Variable types: 5012 continuous, 37 integer (37 binary)

Root relaxation: objective 1.230321e+01, 3899 iterations, 1.28 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   12.30321    0    2          -   12.30321      -     -    2s
H    0     0                       0.0915893   12.30321      -     -    2s
     0     0   12.30321    0    2    0.09159   12.30321      -     -    2s
     0     2   12.30321    0    2    0.09159   12.30321      -     -    3s
*   75    56              27       0.2864668   12.30321  4195%  41.7    4s
    91    62    8.49366   22    5    0.28647   12.30321  4195%  41.9    5s
*  307    74              19       1.2019176   12.29845   923%  58.1    8s
   359    46   12.28821   17    7    1.20192   12.28821   922%  65.3   10s
*  412    15              18       1.4425937    9.52490   560%  68.1   11s

Explored 435 nodes (37087 simplex iterations) in 12.97 seconds
Thread count was 4 (of 4 available processors)

Solution count 4: 1.44259 1.20192 0.286467 0.0915893
Pool objective bound 1.44259

Optimal solution found (tolerance 1.00e-12)
Best objective 1.442593734178e+00, best bound 1.442593734178e+00, gap 0.0000%
optKnock found a optKnock set of large 2 composed by
PFL
and RPE
The production of succinate after optimization is 1.44
The growth rate after optimization is 0.22
The production of other products such as ethanol, formate, lactate and acetate are 15.2, 0.0, -0.0 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.18
The maximun and minimun production of succinate given the optKnock set is 0.06 and 1.22, respectively
...Performing optKnock analysis...
MILP problem with 8400 constraints 63 integer variables and 8370 continuous variables
Optimize a model with 8400 rows, 8370 columns and 34677 nonzeros
Variable types: 8307 continuous, 63 integer (63 binary)
Coefficient statistics:
  Matrix range     [2e-06, 1e+03]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+03]
  RHS range        [1e-02, 1e+03]
Presolve removed 5079 rows and 3321 columns
Presolve time: 0.26s
Presolved: 3321 rows, 5049 columns, 21621 nonzeros
Variable types: 5012 continuous, 37 integer (37 binary)

Root relaxation: objective 1.230321e+01, 3442 iterations, 1.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   12.30321    0    2          -   12.30321      -     -    2s
H    0     0                       0.0915893   12.30321      -     -    2s
H    0     0                       0.0915893   12.30321      -     -    2s
     0     0   12.30321    0    2    0.09159   12.30321      -     -    2s
     0     2   12.30321    0    2    0.09159   12.30321      -     -    2s
*   73    61              22       0.2411614   12.30321  5002%  39.7    4s
*  122    77              19       1.2019176   12.30321   924%  44.1    4s
   133    84   12.27798   14    9    1.20192   12.30321   924%  47.0    5s

Explored 371 nodes (26895 simplex iterations) in 8.23 seconds
Thread count was 4 (of 4 available processors)

Solution count 4: 1.20192 0.241161 0.0915893 0.0915893
Pool objective bound 1.20192

Optimal solution found (tolerance 1.00e-12)
Best objective 1.201917633440e+00, best bound 1.201917633440e+00, gap 0.0000%
optKnock found a optKnock set of large 2 composed by
PFL
and TKT1
The production of succinate after optimization is 1.20
The growth rate after optimization is 0.22
The production of other products such as ethanol, formate, lactate and acetate are 0.0, 0.0, 15.4 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.19
The maximun and minimun production of succinate given the optKnock set is 0.06 and 1.01, respectively

TROUBLESHOOTING 1: "The algorithm takes a long time to find a solution"

TROUBLESHOOTING 2: "The algorithm finds a set of knockouts too big"

TROUBLESHOOTING 3: "The algorithm found a solution that is not useful for me"

EXAMPLE 2: finding optKnock reactions sets of large 3 for increasing production of succinate

% Set optKnock options

% The exchange of succinate will be the objective of the outer problem

options = struct('targetRxn', 'EX_succ_e', 'numDel', 3);

% We will impose that biomass be at least 50% of the biomass of wild-type

constrOpt = struct('rxnList', {{biomass}}, 'values', 0.5*fbaWT.f, 'sense', 'G');

% We will try to find 10 optKnock sets of a maximun length of 3

nIter = 1;

while nIter < threshold

fprintf('...Performing optKnock analysis...')

if isempty(previousSolutions{1})

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt);

else

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions);

end

% determine succinate production and growth rate after optimizacion

succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e'));

growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns, biomass));

etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e'));

formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e'));

lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e'));

acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e'));

setM1 = optKnockSol.rxnList;

if ~isempty(setM1)

previousSolutions{contPreviousSolutions} = setM1;

contPreviousSolutions=contPreviousSolutions + 1;

%printing results

fprintf('optKnock found a optKnock set of large %d composed by ',length(setM1));

for j = 1:length(setM1)

if j == 1

fprintf('%s',setM1{j});

elseif j == length(setM1)

fprintf(' and %s',setM1{j});

else

fprintf(', %s',setM1{j});

end

end

fprintf('\n');

fprintf('The production of succinate after optimization is %.2f \n', succFluxM1);

fprintf('The growth rate after optimization is %.2f \n', growthRateM1);

fprintf('The production of other products such as ethanol, formate, lactate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, lactFluxM1, acetFluxM1);

fprintf('...Performing coupling analysis...\n');

[type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_succ_e');

fprintf('The solution is of type: %s\n', type);

fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth);

fprintf('The maximun and minimun production of succinate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd);

if strcmp(type,'growth coupled')

singleProductionEnvelope(model, setM1, 'EX_succ_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['succ_ex2_' num2str(nIter)], 'outputFolder', 'OptKnockResults');

end

else

if nIter == 1

fprintf('optKnock was not able to found an optKnock set\n');

else

fprintf('optKnock was not able to found additional optKnock sets\n');

end

break;

end

nIter = nIter + 1;

end

...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
ALCD2x
, LDH_D
and PFL
The production of succinate after optimization is 8.81
The growth rate after optimization is 0.15
The production of other products such as ethanol, formate, lactate and acetate are 0.0, 0.0, -0.0 and 4.1, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled
The maximun growth rate given the optKnock set is 0.11
The maximun and minimun production of succinate given the optKnock set is 9.11 and 9.11, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
ACALD
, LDH_D
and PFL
The production of succinate after optimization is 8.81
The growth rate after optimization is 0.15
The production of other products such as ethanol, formate, lactate and acetate are 0.0, 0.0, -0.0 and 4.1, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.11
The maximun and minimun production of succinate given the optKnock set is 4.23 and 9.11, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
LDH_D
, PFL
and TALA
The production of succinate after optimization is 1.20
The growth rate after optimization is 0.22
The production of other products such as ethanol, formate, lactate and acetate are 15.4, 0.0, -0.0 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.19
The maximun and minimun production of succinate given the optKnock set is 0.06 and 1.01, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
GLCabcpp
, PFL
and TALA
The production of succinate after optimization is 1.20
The growth rate after optimization is 0.22
The production of other products such as ethanol, formate, lactate and acetate are 0.0, 0.0, 15.4 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.19
The maximun and minimun production of succinate given the optKnock set is 0.06 and 1.01, respectively

II) LACTATE OVERPRODUCTION

EXAMPLE 1: finding optKnock reactions sets of large 3 for increasing production of lactate

fprintf('\n...EXAMPLE 1: Finding optKnock sets of large 3...\n\n')

...EXAMPLE 1: Finding optKnock sets of large 3...

% Set optKnock options

% The exchange of lactate will be the objective of the outer problem

options = struct('targetRxn', 'EX_lac__D_e', 'numDel', 3);

% We will impose that biomass be at least 50% of the biomass of wild-type

constrOpt = struct('rxnList', {{biomass}}, 'values', 0.5*fbaWT.f, 'sense', 'G');

% We will try to find 10 optKnock sets of a maximun length of 6

previousSolutions = cell(100, 1);

contPreviousSolutions = 1;

nIter = 1;

while nIter < threshold

fprintf('...Performing optKnock analysis...')

if isempty(previousSolutions{1})

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt);

else

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions);

end

% determine lactate production and growth rate after optimizacion

lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e'));

growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns, biomass));

etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e'));

formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e'));

succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e'));

acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e'));

setM1 = optKnockSol.rxnList;

if ~isempty(setM1)

previousSolutions{contPreviousSolutions} = setM1;

contPreviousSolutions=contPreviousSolutions + 1;

%printing results

fprintf('optKnock found a optKnock set of large %d composed by ',length(setM1));

for j = 1:length(setM1)

if j == 1

fprintf('%s', setM1{j});

elseif j == length(setM1)

fprintf(' and %s', setM1{j});

else

fprintf(', %s', setM1{j});

end

end

fprintf('\n');

fprintf('The production of lactate after optimization is %.2f \n', lactFluxM1);

fprintf('The growth rate after optimization is %.2f \n', growthRateM1);

fprintf('The production of other products such as ethanol, formate, succinate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, succFluxM1, acetFluxM1);

fprintf('...Performing coupling analysis...\n');

[type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_lac__D_e');

fprintf('The solution is of type: %s\n', type);

fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth);

fprintf('The maximun and minimun production of lactate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd);

singleProductionEnvelope(model, setM1, 'EX_lac__D_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['lact_ex1_' num2str(nIter)], 'outputFolder', 'OptKnockResults');

else

if nIter == 1

fprintf('optKnock was not able to found an optKnock set\n');

else

fprintf('optKnock was not able to found additional optKnock sets\n');

end

break;

end

nIter = nIter + 1;

end

...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
GLCabcpp
, PFL
and PGI
The production of lactate after optimization is 18.13
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.0, 0.0 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.08
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.72, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
FBA
, PFL
and PGI
The production of lactate after optimization is 18.13
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.0, 0.0 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.08
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.72, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 2 composed by
PFL
and PGI
The production of lactate after optimization is 18.13
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.0, 0.0 and 0.0, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.08
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.72, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 3 composed by
ACKr
, GND
and PGI
The production of lactate after optimization is 18.01
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.6, 0.0 and 0.1, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.09
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.64, respectively

EXAMPLE 2: finding optKnock reactions sets of large 6 for increasing production of lactate

fprintf('...EXAMPLE 3: Finding optKnock sets of large 6...\n')

...EXAMPLE 3: Finding optKnock sets of large 6...

% Set optKnock options

% The exchange of lactate will be the objective of the outer problem

options = struct('targetRxn', 'EX_lac__D_e', 'numDel', 6);

% We will impose that biomass be at least 50% of the biomass of wild-type

constrOpt = struct('rxnList', {{biomass}}, 'values', 0.5*fbaWT.f, 'sense', 'G');

% We will try to find 10 optKnock sets of a maximun length of 2

nIter = 1;

while nIter < threshold

fprintf('...Performing optKnock analysis...')

if isempty(previousSolutions{1})

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt);

else

optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions);

end

% determine lactate production and growth rate after optimizacion

lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e'));

growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns,biomass));

etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e'));

formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e'));

succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e'));

acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e'));

setM1 = optKnockSol.rxnList;

if ~isempty(setM1)

previousSolutions{contPreviousSolutions} = setM1;

contPreviousSolutions = contPreviousSolutions + 1;

%printing results

fprintf('optKnock found a optKnock set of large %d composed by ', length(setM1));

for j = 1:length(setM1)

if j == 1

fprintf('%s', setM1{j});

elseif j == length(setM1)

fprintf(' and %s', setM1{j});

else

fprintf(', %s', setM1{j});

end

end

fprintf('\n');

fprintf('The production of lactate after optimization is %.2f \n', lactFluxM1);

fprintf('The growth rate after optimization is %.2f \n', growthRateM1);

fprintf('The production of other products such as ethanol, formate, succinate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, succFluxM1, acetFluxM1);

fprintf('...Performing coupling analysis...\n');

[type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_lac__D_e');

fprintf('The solution is of type: %s\n', type);

fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth);

fprintf('The maximun and minimun production of lactate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd);

singleProductionEnvelope(model, setM1, 'EX_lac__D_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['lact_ex2_' num2str(nIter)], 'outputFolder', 'OptKnockResults');

else

if nIter == 1

fprintf('optKnock was not able to found an optKnock set\n');

else

fprintf('optKnock was not able to found additional optKnock sets\n');

end

break;

end

nIter = nIter + 1;

end

...Performing optKnock analysis...
optKnock found a optKnock set of large 4 composed by
ACKr
, MDH2, PGI
and RPE
The production of lactate after optimization is 17.96
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.9, 0.0 and 0.1, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.09
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.60, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 5 composed by
ACKr
, GLCabcpp, PGI, RPE
and TALA
The production of lactate after optimization is 17.96
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.9, 0.0 and 0.1, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.09
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.60, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 6 composed by
ALCD2x
, MDH2, MDH3, PGI, PTAr
and RPE
The production of lactate after optimization is 17.96
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.9, 0.0 and 0.1, respectively.
...Performing coupling analysis...
The solution is of type: growth coupled non unique
The maximun growth rate given the optKnock set is 0.09
The maximun and minimun production of lactate given the optKnock set is 10.07 and 18.60, respectively
...Performing optKnock analysis...
optKnock found a optKnock set of large 5 composed by
ACKr
, GLCabcpp, MDH3, PGI
and RPE
The production of lactate after optimization is 17.96
The growth rate after optimization is 0.12
The production of other products such as ethanol, formate, succinate and acetate are 0.0, 0.9, 0.0 and 0.1, respectively.
...Performing coupling analysis...
The solution is of type: non unique
The maximun growth rate given the optKnock set is 0.09
The maximun and minimun production of lactate given the optKnock set is 0.00 and 18.60, respectively

cd(currectDirectory);

TIMING

  1. Example 1 ~ 1-2 minutes
  2. Example 2 ~ 1-2 minutes
  3. Example 3 ~ 1-2 minutes
  4. Example 3 ~ 1-2 minutes
  5. Example 5 ~ 1-2 minutes
  6. Example 6 ~ 1-2 minutes

TROUBLESHOOTING

1) If the algorithm takes a long time to find a solution, it is possible that the seach space is too long. You can reduce the seach space using a smaller set of reactions in the input variable "selectedRxnList"

2) The default number of deletions used by optKnock is 5. If the algorithm is returning more deletions than what you want, you can change the number of deletions using the input variable "numDel"

3) optKnock could find a solution that it is not useful for you. For example, you may think that a solution is very obvious or that it breaks some important biological contraints. If optKnock found a solution that you don't want to find, use the input variable "prevSolutions" to prevent that solution to be found.

ANTICIPATED RESULTS

The optKnock algorithm will find sets of reactions that, when removed from the model, will improve the production of succinate and lactate respectively. In this tutorial, once optKnock finds a solution, then the type of solution is determined (if the product is coupled with biomass formation or not). Some of the sets will generate a coupled solution, i.e., the production rate will increase as biomass formation increases. For these kind of reactions a plot will be generated using the function singleProductionEnvelope and will be saved in the folder tutoriales/optKnock/optKnockResults

When you find a solution with OptKnock, you should always verify the minumum and maximum production rate using the function analizeOptKnock.

References

[1] Burgard, A. P., Pharkya, P. & Maranas, C. D. (2003). OptKnock: A Bilevel Programming Framework for Identifying Gene Knockout Strategies for Microbial Strain Optimization. Biotechnology and Bioengineering, 84(6), 647–657. http://doi.org/10.1002/bit.10803.

[2] Orth, J. D., Conrad, T. M., Na, J., Lerman, J. A., Nam, H., Feist, A. M., & Palsson, B. Ø. (2011). A comprehensive genome‐scale reconstruction of Escherichia coli metabolism—2011. Molecular systems biology, 7(1), 535.

+ \ No newline at end of file diff --git a/tutorials/optKnock/tutorial_optKnock.m b/tutorials/optKnock/tutorial_optKnock.m new file mode 100644 index 0000000000..ef335b123b --- /dev/null +++ b/tutorials/optKnock/tutorial_optKnock.m @@ -0,0 +1,426 @@ +%% OptKnock Tutorial +%% Author: Sebastin N. Mendoza, Center for Mathematical Modeling, University of Chile. snmendoz@uc.cl +%% *Reviewer(s): Thomas Pfau, Sylvian Arreckx, Laurent Heirendt & Ronan Fleming* +%% *INTRODUCTION:* +% In this tutorial we will run optKnock. For a detailed description of the procedure, +% please see [1]. Briefly, the problem is to find a set of reactions of size "numDel" +% such that when these reactions are deleted from the network, the mutant created +% will produce a particular target of interest in a higher rate than the wild-type +% strain. +% +% For example, imagine that we would like to increase the production of succinate +% or lactate in Escherichia coli. Which are the knock-outs needed to increase +% the production of these products? We will approach these problems in this tutorial. +%% MATERIALS +%% EQUIPMENT +% # MATLAB +% # A solver for Mixed Integer Linear Programming (MILP) problems. For example, +% Gurobi. +%% *EQUIPMENT SETUP* +% Use changeCobraSolver to choose the solver for MILP problems. +%% PROCEDURE +% The proceduce consists on the following steps +% +% 1) Define contraints +% +% 2) Define a set of reactions to search knockouts. Only reactions in this +% set will be deleted. +% +% 3) Define the number of reactions to be deleter, the target reaction and +% some constraints to be accomplish +% +% 4) Run optKnock. *TIMING: *This task should take from a few seconds to +% a few hours depending on the size of your reconstruction +% +% We verify that cobratoolbox has been initialized and that the solver has +% been set. + +global TUTORIAL_INIT_CB; +if ~isempty(TUTORIAL_INIT_CB) && TUTORIAL_INIT_CB==1 + initCobraToolbox + changeCobraSolver('gurobi','all'); +end + +fullPath = which('optKnockTutorial'); +folder = fileparts(fullPath); +currectDirectory = pwd; +cd(folder); + +%% +% We load the model of E. coli [2] + +model = readCbModel('iJO1366.mat') +biomass = 'BIOMASS_Ec_iJO1366_core_53p95M'; +%% +% We define the maximum number of solutions to find + +threshold = 5; +%% +% First, we define the set for reactions which could be deleted from the +% network. Reactions not in this list are not going to be deleted. + +selectedRxnList = {'GLCabcpp'; 'GLCptspp'; 'HEX1'; 'PGI'; 'PFK'; 'FBA'; 'TPI'; 'GAPD'; 'PGK'; 'PGM'; 'ENO'; 'PYK'; 'LDH_D'; 'PFL'; 'ALCD2x'; 'PTAr'; 'ACKr'; 'G6PDH2r'; 'PGL'; 'GND'; 'RPI'; 'RPE'; 'TKT1'; 'TALA'; 'TKT2'; 'FUM'; 'FRD2'; 'SUCOAS'; 'AKGDH'; 'ACONTa'; 'ACONTb'; 'ICDHyr'; 'CS'; 'MDH'; 'MDH2'; 'MDH3'; 'ACALD'}; + +%% +% Then, we define some constraints + + +% prespecified amount of glucose uptake 10 mmol/grDW*hr +model = changeRxnBounds(model, 'EX_glc__D_e', -10, 'b'); + +% Unconstrained uptake routes for inorganic phosphate, sulfate and +% ammonia +model = changeRxnBounds(model, 'EX_o2_e', 0, 'l'); +model = changeRxnBounds(model, 'EX_pi_e', -1000, 'l'); +model = changeRxnBounds(model, 'EX_so4_e', -1000, 'l'); +model = changeRxnBounds(model, 'EX_nh4_e', -1000, 'l'); + +% The optimization step could opt for or against the phosphotransferase +% system, glucokinase, or both mechanisms for the uptake of glucose +model = changeRxnBounds(model, 'GLCabcpp', -1000, 'l'); +model = changeRxnBounds(model, 'GLCptspp', -1000, 'l'); +model = changeRxnBounds(model, 'GLCabcpp', 1000, 'u'); +model = changeRxnBounds(model, 'GLCptspp', 1000, 'u'); +model = changeRxnBounds(model, 'GLCt2pp', 0, 'b'); + +% Secretion routes for acetate, carbon dioxide, ethanol, formate, lactate +% and succinate are enabled +model = changeRxnBounds(model, 'EX_ac_e', 1000, 'u'); +model = changeRxnBounds(model, 'EX_co2_e', 1000, 'u'); +model = changeRxnBounds(model, 'EX_etoh_e', 1000, 'u'); +model = changeRxnBounds(model, 'EX_for_e', 1000, 'u'); +model = changeRxnBounds(model, 'EX_lac__D_e', 1000, 'u'); +model = changeRxnBounds(model, 'EX_succ_e', 1000, 'u'); +%% +% Then, we calculates the production of metabolites before running optKnock + +% determine succinate production and growth rate before optimizacion +fbaWT = optimizeCbModel(model); +succFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_succ_e')); +etohFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_etoh_e')); +formFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_for_e')); +lactFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_lac__D_e')); +acetFluxWT = fbaWT.x(strcmp(model.rxns, 'EX_ac_e')); +growthRateWT = fbaWT.f; +fprintf('The production of succinate before optimization is %.1f \n', succFluxWT); +fprintf('The growth rate before optimization is %.1f \n', growthRateWT); +fprintf('The production of other products such as ethanol, formate, lactate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxWT, formFluxWT, lactFluxWT, acetFluxWT); +%% +% *I) SUCCINATE OVERPRODUCTION* +% +% *EXAMPLE 1:* *finding optKnock reactions sets of large 2 for increasing +% production of succinate* + +fprintf('\n...EXAMPLE 1: Finding optKnock sets of large 2 or less...\n\n') +% Set optKnock options +% The exchange of succinate will be the objective of the outer problem +options = struct('targetRxn', 'EX_succ_e', 'numDel', 2); +% We will impose that biomass be at least 50% of the biomass of wild-type +constrOpt = struct('rxnList', {{biomass}},'values', 0.5*fbaWT.f, 'sense', 'G'); +% We will try to find 10 optKnock sets of a maximun length of 2 +previousSolutions = cell(10, 1); +contPreviousSolutions = 1; +nIter = 1; +while nIter < threshold + fprintf('...Performing optKnock analysis...\n') + if isempty(previousSolutions{1}) + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt); + else + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions, 1); + end + + % determine succinate production and growth rate after optimizacion + succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e')); + growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns, biomass)); + etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e')); + formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e')); + lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e')); + acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e')); + setM1 = optKnockSol.rxnList; + + if ~isempty(setM1) + previousSolutions{contPreviousSolutions} = setM1; + contPreviousSolutions = contPreviousSolutions + 1; + %printing results + fprintf('optKnock found a optKnock set of large %d composed by ', length(setM1)); + for j = 1:length(setM1) + if j == 1 + fprintf('%s', setM1{j}); + elseif j == length(setM1) + fprintf(' and %s', setM1{j}); + else + fprintf(', %s', setM1{j}); + end + end + fprintf('\n'); + fprintf('The production of succinate after optimization is %.2f \n', succFluxM1); + fprintf('The growth rate after optimization is %.2f \n', growthRateM1); + fprintf('The production of other products such as ethanol, formate, lactate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, lactFluxM1, acetFluxM1); + fprintf('...Performing coupling analysis...\n'); + [type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_succ_e'); + fprintf('The solution is of type: %s\n', type); + fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth); + fprintf('The maximun and minimun production of succinate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd); + if strcmp(type, 'growth coupled') + singleProductionEnvelope(model, setM1, 'EX_succ_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['succ_ex1_' num2str(nIter)], 'outputFolder', 'OptKnockResults'); + end, + else + if nIter == 1 + fprintf('optKnock was not able to found an optKnock set\n'); + else + fprintf('optKnock was not able to found additional optKnock sets\n'); + end + break; + end + nIter = nIter + 1; +end + + +%% +% *TROUBLESHOOTING 1*: "The algorithm takes a long time to find a solution" +% +% *TROUBLESHOOTING 2*: "The algorithm finds a set of knockouts too big" +% +% *TROUBLESHOOTING 3*: "The algorithm found a solution that is not useful +% for me" +% +% *EXAMPLE 2:* *finding optKnock reactions sets of large 3 for increasing +% production of succinate* + + +% Set optKnock options +% The exchange of succinate will be the objective of the outer problem +options = struct('targetRxn', 'EX_succ_e', 'numDel', 3); +% We will impose that biomass be at least 50% of the biomass of wild-type +constrOpt = struct('rxnList', {{biomass}}, 'values', 0.5*fbaWT.f, 'sense', 'G'); +% We will try to find 10 optKnock sets of a maximun length of 3 +nIter = 1; +while nIter < threshold + fprintf('...Performing optKnock analysis...') + if isempty(previousSolutions{1}) + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt); + else + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions); + end + + % determine succinate production and growth rate after optimizacion + succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e')); + growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns, biomass)); + etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e')); + formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e')); + lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e')); + acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e')); + setM1 = optKnockSol.rxnList; + + if ~isempty(setM1) + previousSolutions{contPreviousSolutions} = setM1; + contPreviousSolutions=contPreviousSolutions + 1; + %printing results + fprintf('optKnock found a optKnock set of large %d composed by ',length(setM1)); + for j = 1:length(setM1) + if j == 1 + fprintf('%s',setM1{j}); + elseif j == length(setM1) + fprintf(' and %s',setM1{j}); + else + fprintf(', %s',setM1{j}); + end + end + fprintf('\n'); + fprintf('The production of succinate after optimization is %.2f \n', succFluxM1); + fprintf('The growth rate after optimization is %.2f \n', growthRateM1); + fprintf('The production of other products such as ethanol, formate, lactate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, lactFluxM1, acetFluxM1); + fprintf('...Performing coupling analysis...\n'); + [type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_succ_e'); + fprintf('The solution is of type: %s\n', type); + fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth); + fprintf('The maximun and minimun production of succinate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd); + if strcmp(type,'growth coupled') + singleProductionEnvelope(model, setM1, 'EX_succ_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['succ_ex2_' num2str(nIter)], 'outputFolder', 'OptKnockResults'); + end + else + if nIter == 1 + fprintf('optKnock was not able to found an optKnock set\n'); + else + fprintf('optKnock was not able to found additional optKnock sets\n'); + end + break; + end + nIter = nIter + 1; +end + +%% +% *II) LACTATE OVERPRODUCTION* +% +% *EXAMPLE 1: finding optKnock reactions sets of large 3 for increasing production +% of lactate* + +fprintf('\n...EXAMPLE 1: Finding optKnock sets of large 3...\n\n') +% Set optKnock options +% The exchange of lactate will be the objective of the outer problem +options = struct('targetRxn', 'EX_lac__D_e', 'numDel', 3); +% We will impose that biomass be at least 50% of the biomass of wild-type +constrOpt = struct('rxnList', {{biomass}}, 'values', 0.5*fbaWT.f, 'sense', 'G'); +% We will try to find 10 optKnock sets of a maximun length of 6 +previousSolutions = cell(100, 1); +contPreviousSolutions = 1; +nIter = 1; +while nIter < threshold + fprintf('...Performing optKnock analysis...') + if isempty(previousSolutions{1}) + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt); + else + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions); + end + + % determine lactate production and growth rate after optimizacion + lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e')); + growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns, biomass)); + etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e')); + formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e')); + succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e')); + acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e')); + setM1 = optKnockSol.rxnList; + + if ~isempty(setM1) + previousSolutions{contPreviousSolutions} = setM1; + contPreviousSolutions=contPreviousSolutions + 1; + %printing results + fprintf('optKnock found a optKnock set of large %d composed by ',length(setM1)); + for j = 1:length(setM1) + if j == 1 + fprintf('%s', setM1{j}); + elseif j == length(setM1) + fprintf(' and %s', setM1{j}); + else + fprintf(', %s', setM1{j}); + end + end + fprintf('\n'); + fprintf('The production of lactate after optimization is %.2f \n', lactFluxM1); + fprintf('The growth rate after optimization is %.2f \n', growthRateM1); + fprintf('The production of other products such as ethanol, formate, succinate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, succFluxM1, acetFluxM1); + fprintf('...Performing coupling analysis...\n'); + [type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_lac__D_e'); + fprintf('The solution is of type: %s\n', type); + fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth); + fprintf('The maximun and minimun production of lactate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd); + singleProductionEnvelope(model, setM1, 'EX_lac__D_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['lact_ex1_' num2str(nIter)], 'outputFolder', 'OptKnockResults'); + else + if nIter == 1 + fprintf('optKnock was not able to found an optKnock set\n'); + else + fprintf('optKnock was not able to found additional optKnock sets\n'); + end + break; + end + nIter = nIter + 1; +end + +%% +% +% +% *EXAMPLE 2: finding optKnock reactions sets of large 6 for increasing production +% of lactate* + +fprintf('...EXAMPLE 3: Finding optKnock sets of large 6...\n') +% Set optKnock options +% The exchange of lactate will be the objective of the outer problem +options = struct('targetRxn', 'EX_lac__D_e', 'numDel', 6); +% We will impose that biomass be at least 50% of the biomass of wild-type +constrOpt = struct('rxnList', {{biomass}}, 'values', 0.5*fbaWT.f, 'sense', 'G'); +% We will try to find 10 optKnock sets of a maximun length of 2 +nIter = 1; +while nIter < threshold + fprintf('...Performing optKnock analysis...') + if isempty(previousSolutions{1}) + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt); + else + optKnockSol = OptKnock(model, selectedRxnList, options, constrOpt, previousSolutions); + end + + % determine lactate production and growth rate after optimizacion + lactFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_lac__D_e')); + growthRateM1 = optKnockSol.fluxes(strcmp(model.rxns,biomass)); + etohFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_etoh_e')); + formFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_for_e')); + succFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_succ_e')); + acetFluxM1 = optKnockSol.fluxes(strcmp(model.rxns, 'EX_ac_e')); + setM1 = optKnockSol.rxnList; + + if ~isempty(setM1) + previousSolutions{contPreviousSolutions} = setM1; + contPreviousSolutions = contPreviousSolutions + 1; + %printing results + fprintf('optKnock found a optKnock set of large %d composed by ', length(setM1)); + for j = 1:length(setM1) + if j == 1 + fprintf('%s', setM1{j}); + elseif j == length(setM1) + fprintf(' and %s', setM1{j}); + else + fprintf(', %s', setM1{j}); + end + end + fprintf('\n'); + fprintf('The production of lactate after optimization is %.2f \n', lactFluxM1); + fprintf('The growth rate after optimization is %.2f \n', growthRateM1); + fprintf('The production of other products such as ethanol, formate, succinate and acetate are %.1f, %.1f, %.1f and %.1f, respectively. \n', etohFluxM1, formFluxM1, succFluxM1, acetFluxM1); + fprintf('...Performing coupling analysis...\n'); + [type, maxGrowth, maxProd, minProd] = analyzeOptKnock(model, setM1, 'EX_lac__D_e'); + fprintf('The solution is of type: %s\n', type); + fprintf('The maximun growth rate given the optKnock set is %.2f\n', maxGrowth); + fprintf('The maximun and minimun production of lactate given the optKnock set is %.2f and %.2f, respectively \n\n', minProd, maxProd); + singleProductionEnvelope(model, setM1, 'EX_lac__D_e', biomass, 'savePlot', 1, 'showPlot', 1, 'fileName', ['lact_ex2_' num2str(nIter)], 'outputFolder', 'OptKnockResults'); + else + if nIter == 1 + fprintf('optKnock was not able to found an optKnock set\n'); + else + fprintf('optKnock was not able to found additional optKnock sets\n'); + end + break; + end + nIter = nIter + 1; +end + +cd(currectDirectory); +%% TIMING +% # Example 1 ~ 1-2 minutes +% # Example 2 ~ 1-2 minutes +% # Example 3 ~ 1-2 minutes +% # Example 3 ~ 1-2 minutes +% # Example 5 ~ 1-2 minutes +% # Example 6 ~ 1-2 minutes +%% TROUBLESHOOTING +% 1) If the algorithm takes a long time to find a solution, it is possible that +% the seach space is too long. You can reduce the seach space using a smaller +% set of reactions in the input variable "selectedRxnList" +% +% 2) The default number of deletions used by optKnock is 5. If the algorithm +% is returning more deletions than what you want, you can change the number of +% deletions using the input variable "numDel" +% +% 3) optKnock could find a solution that it is not useful for you. For example, +% you may think that a solution is very obvious or that it breaks some important +% biological contraints. If optKnock found a solution that you don't want to find, +% use the input variable "prevSolutions" to prevent that solution to be found. +%% ANTICIPATED RESULTS +% The optKnock algorithm will find sets of reactions that, when removed from +% the model, will improve the production of succinate and lactate respectively. +% In this tutorial, once optKnock finds a solution, then the type of solution +% is determined (if the product is coupled with biomass formation or not). Some +% of the sets will generate a coupled solution, i.e., the production rate will +% increase as biomass formation increases. For these kind of reactions a plot +% will be generated using the function singleProductionEnvelope and will be saved +% in the folder tutoriales/optKnock/optKnockResults +% +% When you find a solution with OptKnock, you should always verify the minumum +% and maximum production rate using the function analizeOptKnock. +%% References +% [1] Burgard, A. P., Pharkya, P. & Maranas, C. D. (2003). OptKnock: A Bilevel +% Programming Framework for Identifying Gene Knockout Strategies for Microbial +% Strain Optimization. Biotechnology and Bioengineering, 84(6), 647657. http://doi.org/10.1002/bit.10803. +% +% [2] Orth, J. D., Conrad, T. M., Na, J., Lerman, J. A., Nam, H., Feist, +% A. M., & Palsson, B. . (2011). A comprehensive genome?scale reconstruction +% of Escherichia coli metabolism2011. _Molecular systems biology_, _7_(1), 535. \ No newline at end of file