BACKGROUND
metabolic engineering aims to improve product yields in cellular systems by applying a variety of tools.
constraint based methods which use only the stoichiometry of metabolic reactions have been particularly successful in the development of strategies towards fulfilling this aim  <cit> .
one important application is the prediction of knockouts to enforce desired metabolic behaviors in an organism.
a method that allows one to predict efficient intervention strategies using the concept of minimal cut sets mcss, was developed by klamt and gilles  <cit> .
this was generalized to constrained minimal cut sets cmcs, where in addition to blocking undesired fluxes, survival of some desired fluxes is possible  <cit> .
the automatic partitioning method apm uses an objective function to specify the design objectives and the partitioning of fluxes into desired/undesired is done automatically to find successively larger cmcs till a global optimum is reached  <cit> .
previously we showed that a genetic algorithm could reach the global optimum faster than than apm  <cit> .
however, all these methods are applicable only to small and medium-scale metabolic networks.
in a recent work by ballerstein et al., it was shown that cmcs can be directly calculated from the stoichiometric matrix  <cit> .
using this method, it is possible to calculate intervention strategies even in genome-scale metabolic networks  <cit> .
another work extended this concept to include regulation  <cit> .
a limitation of this method is that the desired flux or flux ratio of a metabolite has to be manually specified to get corresponding cmcs.
there exist other constraint based methods for predicting intervention strategies.
optknock solves a bi-level optimization problem, to predict knockouts leading to maximal product formation at maximal growth  <cit> .
a three-level optimization problem is used to maximize minimal product formation in robustknock  <cit> .
optgene uses a genetic algorithm to predict knockouts  <cit> .
similarly, evolutionary algorithms and simulated annealing have been used in  <cit> .
another metaheuristic approach was using a hybrid of bees algorithm with flux balance analysis fba  <cit> .
while these methods optimize for design goals, doing so with a minimal number of knockouts is not necessarily guaranteed.
from an engineering perspective, we would like the organism to have a guaranteed high yield for the product of interest.
given that even in the face of genetic perturbations microorganisms redirect metabolic flux towards maximizing cellular growth  <cit> , this high yield must be maintained at high growth rates.
additionally, the number of knockouts should be as small as possible to facilitate easy implementation in the laboratory.
here we present a new method, psomcs, which uses particle swarm optimization pso along with the method developed in  to calculate cmcs while overcoming the mentioned limitations of other methods.
our basic motivation is to combine the computational rigour of cmcs with the flexibility of the optimization-based approaches in order to solve  intervention problems efficiently.
we aim to find not only the optimal intervention strategy for a given design but also the best possible design.
in addition, we show that psomcs is also faster than other methods which try to find cmcs leading to optimal design objectives.
methods
calculating cmcs
a metabolic network of m internal metabolites connected by n reactions in steady state is represented by the set of linear equations
<dig> nr= <dig>
where n is a m×n matrix consisting of stoichiometric coefficients of all participating reactions such that each column represents one reaction.
r is a vector of reaction fluxes.
reactions can be both reversible  and irreversible , thereby imposing the constraint
<dig> ri≥0∀i∈irrev.
and  define a flux space.
depending on the desired outcome, an intervention problem can be set up dividing this space into desired and undesired fluxes.
the set of undesired fluxes for t reactions can be defined by
<dig> tr≤t
where t∈ℝt×n and t∈ℝt× <dig>  likewise, the set of desired fluxes for d reactions can be defined by
<dig> dr≤d
with d∈ℝd×n and d∈ℝd× <dig>
in  <cit> , cmcs are calculated by first solving a series of mixed integer linear programming milp problems representing  and  and then filtering those solutions which also satisfy .
in  <cit> , this is combined into a single system represented as  in  <cit> )
<dig> nrevtirev−irevtrevt0nirrtiirr−iirrtirrt00000n0000d×uvpvnwr=≥=≤000dttw≤−cu∈ℝm,vp,vn∈ℝn,d∈ℝd,vp,vn,w,rirr≥ <dig> c> <dig>
note that the n and t matrices have been split into reversible  and irreversible submatrices .
similarly, identity submatrices for reversible and irreversible reactions are represented by the matrices i
rev and i
irr respectively.
cmcs are directly calculated by finding solutions with minimum number of non-zero entries in v
p,v
n. additionally binary indicator variables z
p and z
n are introduced such that z
p
i= <dig> if v
p
i= <dig> and z
p
i= <dig> if v
p
i> <dig> and similarly for zn, vn.
only one direction of v  can be active, hence
<dig> zpi+zni≤ <dig>
we set up the following optimization problem
<dig> minimize∑i=1ns.t.,
with the additional constraint that the flux through a reaction is turned off if it is part of a cmcs, i.e., r
i= <dig> if z
p
i= <dig> || z
n
i= <dig>
with this system it is possible to find cmcs which will result in designs satisfying constraints on yields/fluxes specified by , .
however, we would like to have a method which given some design objectives  calculates cmcs corresponding to optimal values for the design objectives.
since any design can be represented as a function of t,d,t and d, the optimization problem can be stated as
<dig> maxfs.t..
in other words, the problem is to find optimal combinations of {target/desired} yields for all reactions to be optimized.
this is not easy for a few reasons.
in general, this is a non-linear optimization problem.
non-linear optimization is known to be inherently complex with general deterministic solutions being impossible to find.
secondly, slight adjustments in ,  could result in completely different cmcs with different cardinalities.
finally, not all such combinations will result in cmcs.
these issues become acute when the search space is more dense with many possible combinations, as in large and genome-scale metabolic networks.
we attack this problem using pso as it has been successfully used to find solutions to complex non-linear optimization problems in other fields .
particle swarm optimization
pso is a metaheuristic inspired by the flocking behavior of birds  <cit> .
in pso, particles distributed within a multidimensional space collectively move towards an optimum guided by a fitness function.
particle fitness is determined by its position in the search space.
the motion of a particle is influenced by its neighbours and the currently known fittest particle.
more information on pso can be found in .
typically, a particle is made up of three j-dimensional vectors, where j is the dimensionality of the search space.
these represent the current position x, its previous best position p which is the position corresponding to the highest fitness achieved by the particle and the velocity v, fig.
<dig>  particle motion is guided by the following equations,
<dig> vi=χ{vi+φ1β1+φ2β2}
fig.
<dig> schematic of the pso particle.
a particle stores three types of information: the current values, values corresponding to its own previous best fitness and velocities corresponding to each objective
<dig> xi=xi+vi
i∈{ <dig> .j}.
g is the position corresponding to the global best fitness of the entire swarm till the current t. φ
<dig>  are called “acceleration constants” and determine the relative influence of the particle’s own knowledge and that of the group, both of which are commonly set to  <dig>  <cit> .
β
<dig>  are uniformly generated random numbers within the range  needs constraints specified by ,  to calculate corresponding cmcs.
for example, consider a network which has, among other reactions, a substrate uptake reaction r
s, a reaction for the product secretion r
p and one for biomass r
bio.
an optimal design could be stated as having r
p/r
s≥x
<dig> and also that biomass fluxes of r
bio/r
s≥x
<dig> exist.
however, we don’t know the combinations of x
<dig>  x
<dig> resulting in optimal design.
this is where a pso can be useful.
after initializing x
,
v, the set of positions and velocities for all particles, within the range of values for {x
<dig> x
2} on some constant r
s, the pso iteratively finds increasingly better solutions for  using  and  and moves towards the global optimum.
the psomcs flowchart is shown in fig.
<dig>
fig.
<dig> flowchart of psomcs.
p and g are the current particle best and global best respectively.
the algorithm stops when the number of iterations reaches a pre-specified maximum or if the maximum fitness remains unchanged for a pre-specified number of iterations
the fitness function will depend on the nature of the desired optimum.
considering that our objective is to have a design with high yields and minimal knockouts, the following fitness function was used,
<dig> f=1−|cmcs|n·∏ixixi.
RESULTS
to clarify the working of psomcs, we first apply our method to a small toy network, optimizing for only a single reaction.
next, to confirm the accuracy of our predictions, we compare our method against another method based on a genetic algorithm  which we had previously developed  <cit> .
the model used is the medium-scale e. coli core model presented in  <cit> .
finally we find optimal intervention strategies for maximizing the minimal product yield in a genome-scale metabolic network.
fba was used to calculate the range of yields  for each objective and particles were initialised within this range.
only one solution is calculated for a milp.
the parameters used are shown in table  <dig>  implementation of psomcs was done using perl http://www.perl.org/.
for the performance critical parts of the program, i.e., solving the milp and also the lp, the ibm ilog cplex optimization studio - a commercial optimization package - was used through the math::cplex perl module.
also, our algorithm is designed to make use of modern cpu architectures and can be run in parallel on multiple cores.
details of parameters used for the different models
consider the network given in fig.
<dig>  we wish to find minimal knockouts which will ensure the highest possible yield for reaction r <dig>  in the first iteration, cmcs corresponding to low yields are found.
in the second iteration, all particles move towards higher yields.
one particle, on the solution of its dual system gives the cmcs of ‘r <dig> r9’.
removal of r <dig> and r <dig> from the network blocks all flux through r <dig> and r <dig>  thus redirecting the network flux through r <dig>  this corresponds to the highest minimal yield of  <dig> for r <dig>
fig.
<dig> psomcs small example.
running the psomcs on a toy network.
this network has three input reactions, which can be assumed to be substrates and three secretion reactions, which can be assumed to be three different products.
we want to maximise the yield of r <dig>  that is maximize ).
note that the particles operate in a single dimensional search space and x represents the yield for r <dig>  after performing fba to determine the maximum and minimum yields for r <dig> given unit substrate uptake, four particles are initialised within this range.
initial velocities are also assigned.
cmcss are calculated after creating and solving the dual system.
fitness is a function of x and the cardinality of the cmcs.
g corresponds to x with the highest fitness which is particle  <dig> after both the first and second iterations.
after the first iteration, every particle except the first has a value for p. note that for particle  <dig> a yield higher than  <dig>  is guaranteed.
in reality, the minimal yield with the corresponding cmcs is  <dig>  which is also the case for particle  <dig>  this is the value the algorithm will return if allowed to run for a few more iterations
we apply psomcs to generate designs in an e. coli core network which will ensure high yield of ethanol even in the face of high growth.
this network was previously used to design a high yield ethanol producing strain in  <cit> .
this model has  <dig> reactions and  <dig> metabolites.
we had previously used this model to predict optimal intervention strategies using a genetic algorithm , which we had shown to be faster than other current approaches  <cit> , particularly compared to apm, which is guaranteed to find the optimal solution  <cit> .
here we compare our approach with gamcs in terms of speed and accuracy of results.
the machine used had the following specifications −  <dig> cpus,  <dig> cores, intel xeon x <dig>  <dig>  ghz, running an ubuntu  <dig>  lts operating system.
the time taken for a typical psomcs and gamcs run is plotted in fig.
4
a. the superiority of our method in terms of speed can be clearly observed.
gamcs takes  <dig>  seconds to reach the maximum fitness.
psomcs takes only  <dig> seconds for the same.
this is an over  <dig> fold improvement in performance.
in comparison, apm would not only require that the desired efms be assigned weights, but also the time taken by it would have been outside the boundaries of this plot.
the cmcs corresponding to the optimum obtained by both gamcs and psomcs are exactly the same.
figure 4
b is one of the designs corresponding to a high fitness.
this design was in the solution pool of both the pso and ga methods.
in this design, a minimum ethanol yield of  <dig>  is guaranteed even when the growth rate is  <dig> .
also, as can be expected, production of competing by-products: acetate, lactate and succinate is blocked.
additionally, flux through the oxidative part of the pentose phosphate pathway is blocked and so is the pyruvate-malate cycling.
multiple cmcs resulting in similar design characteristics were returned by our method.
fig.
<dig> comparing the runtimes of psomcs and gamcs.
a plotting the runtimes of psomcs and the ga we had previously implemented clearly shows psomcs is orders of magnitude faster than gamcs.
note that the time axis is logarithmic and that both algorithms reach the same maximum fitness.
the change in knockout sizes are indecated by the numbers along the line.
b both methods also produce similar designs, an example of which is shown.
this design is obtained with  <dig> knockouts .
the plot was generated by applying the knockouts on the complete set of  <dig> efms of the escherichia coli core model.
r_norm is the sum of uptake rates for the five carbon substrates, glucose, galactose, mannose, arabinose and xylose under aerobic conditions
to test the capabilities of our method we applied it to the genome-scale model of e. coli presented in  <cit> .
our aim was to find cmcs that result in an scenario of growth-coupled ethanol yield.
a few strategies were used in  <cit>  to reduce the network size.
these strategies are aimed at reducing the network size and improving computational efficiency, which takes real growth conditions into account and removing all superfluous components.
first, the network was reduced to grow anaerobically on glucose as the only carbon source.
the resulting network has  <dig> reactions and  <dig> metabolites.
network compression was done by combining reactions operating at fixed ratios into reaction subsets.
exchange reactions, spontaneous reactions and reactions essential for the ethanol and biomass production were excluded from participating in cmcs by setting their corresponding zp, zn variables to zero.
the machine we used for this test had  <dig> cpus, 396gb ram, intel xeon e5- <dig>  <dig>  ghz processor, running on ubuntu  <dig>  lts.
the cmcs cardinality was limited to  <dig>  with  <dig> particles being processed in parallel, the program was run for  <dig> iterations.
it took  <dig> iterations  to find the optimal design.
one of the designs is shown in fig.
<dig> along with designs obtained using optknock and robustknock on the same machine.
the envelope of the strain specific phenotypic solution space was calculated with flux variability analysis fva  <cit>  of the iaf <dig> network while considering the respective knockouts predicted by each method.
the minimally required biomass production was set at  <dig>  and both were limited by unit glucose uptake and a maximum knockout size of  <dig>  optknock took  <dig> minutes to run while robustknock ran for  <dig> minutes.
the minimal ethanol yields were  <dig> in both cases.
as can be observed, psomcs offers a better design with the ethanol production being strongly coupled to biomass production and at no point falls below a yield of  <dig> .
fig.
<dig> design for a genome-scale e. coli model.
e. coli was designed for enhanced ethanol production using the genome-scale iaf <dig> model.
for comparison, designs obtained using optknock and robustknock are also presented.
the design using psomcs guarantees a minimal ethanol yield of  <dig> , in contrast this is  <dig> for both robustknock and optknock.
all designs have a maximum biomass production rate greater than  <dig>  with the one for psomcs being comparatively lower.
the maximum yield for all the designs is  <dig>  the given plots have been generated by using fva on the iaf <dig> model while considering the respective knockouts produced by each method.
the fba solution space at maximum growth is highlighted, with crosses indicating the maximum and squares the minimum ethanol yield.
all designs involve  <dig> knockouts -  for psomcs,  for optknock and  for robustknock
discussion
here we have presented a method, psomcs, to design strains with high minimal product yield using knockouts of minimal possible size.
to do this, we employ a pso together with the direct enumeration of cmcs developed in .
this method has made it possible to find cmcs in large and genome-scale networks.
however, it is not designed to optimize engineering goals.
that is, we we would like to find not only the optimal intervention strategy for a given design but the best possible design too.
finding intervention strategies that achieve this is an important goal of metabolic engineering, especially in the production of industrially important chemicals.
we deliver on this goal by using a pso built on top of the base provided by the direct enumeration of cmcs.
our method thus expands the utility of this method.
additionally we would like to point out that in the case of optimizing for a single reaction, solving  with continuous values within the  range for that reaction would suffice.
however, in the presence of multiple objectives this task becomes computationally exhaustive and infeasible, thereby justifying the use of a metaheuristic approach such as the one used here.
there have been other methods with a similar strategy as ours, which is the use of a metaheuristic in combination with another method like linear programming.
most methods have relied on genetic algorithms  <cit> , evolutionary algorithms and simulated annealing  <cit>  and also an artificial bees algorithm  <cit> .
ours is the first attempt at using the dual method in a similar fashion, along with the use of a pso.
as shown by the comparison with optknock and robustknock in fig.
<dig>  although all designs have the same highest ethanol yield of  <dig>  psomcs provides a design with the highest guaranteed minimal ethanol yield.
robustknock was developed to overcome the ’too-optimistic’ nature of optknock and this is reflected in the nature of their respective designs.
also of note is the fact that both optknock and robustknock need a minimal level of biomass production to be manually specified while psomcs does not.
in fact, if we reduce the minimal biomass production requirement to  <dig>  , robustknock runs for over  <dig> hours without finding the optimum.
running optknock and robustknock multiple times with different biomass levels will result in different solutions, some of which will be better than others.
psomcs eliminates this need to manually set reaction fluxes and searches the entire feasible space of biomass yields to find the optimal one.
growth-coupling is a key principle in metabolic engineering.
it requires that growth should only be feasible if a desired compound, like ethanol, is mandatorily produced as by-product.
it can be seen in fig.
<dig> that psomcs achieves this with a growth rate about one third of the wild-type.
however, growth-coupling does not enforce nor require that the maximal product yield is attained at a non-zero growth rate.
in fact fig.
<dig> illustrates the rule rather than the exception, as typically the maximum product yield is achieved at zero growth  <cit> .
furthermore, an ideal production state will be characterized by zero growth, where all available resources are used for product formation.
in this senses, biomass production can be seen as an “unwanted” by-product.
recent advances in fermentation processes employ zero-growth approaches  <cit> .
however, these approaches are associated with many challenges which go far beyond the scope of the presented work.
nevertheless, fig.
<dig> indicates that the presented designs retain their wild-type behavior to be operated as optimal zero-growth factories.
in heuristic search algorithms, performance comes at the cost of being too specific to the problem being solved  <cit> .
by virtue of having few parameters, psos are less affected by this problem.
in our implementation, we have used parameter values as found in the general pso literature without the need to adjust them.
the only parameters that we adjusted were the number of particles and the number of iterations.
we clearly use fewer particles than is typical.
this is because we found a population size of  <dig> to be sufficient for our needs .
although we have sampled the entire solution space, particles can easily be forced to explore a subspace.
certain reactions can be excluded from being considered for knockouts by forcing their corresponding indicator variables in the dual system to be  <dig>  our fitness function is specific to our target design, however new fitness functions can be thought of depending on the desired final objective.
our method produces cmcs leading to designs with similar characteristics as the one used in  <cit> .
our method also returns multiple solutions.
the limiting factor in our method is the milp for the dual system.
milps are more difficult to solve than lps and may consume large amounts of time as well as memory  <cit> .
during our runs, the search tree generated by cplex’s branch and cut algorithm for a single milp grew to consume over  <dig> gb of memory when limited to a knockout size of  <dig>  this memory consumption grows quickly with increasing knockout size, thereby limiting the ability of psomcs to find the optimal solution.
improvements in run time can be made by forcing psomcs to explore only a part of the flux space leading to a smaller solution space to be explored.
for instance lets consider the design in fig.
<dig>  with a minimal biomass yield of  <dig> , the optimal design presented here was found within  <dig> hours.
further improvements to performance could be obtained by following the strategies outlined in  <cit> .
also, algorithmic improvements in solving milps could be useful in this regard.
here we have dealt only with knockout strategies to design better strains.
it can easily be extended to include the concept of regulatory mcs introduced in  <cit>  which combine reaction up/downregulation with knockouts.
there are other constraint based methods dealing with intervention strategies like gene knock-ins and up/downregulation.
psos and swarm intelligence algorithms in general may be used to compliment these methods.
CONCLUSIONS
psomcs finds the best possible design in metabolic networks given multiple objectives with the corresponding cmcs.
we have demonstrated its capability in finding optimal knockouts and designs in genome-scale metabolic networks.
it finds competitive designs compared to standard tools and is orders of magnitude faster than efm based tools in finding the optimal solution.
psomcs could be used to predict minimal knockouts resulting in optimal yields in industrially important microorganisms.
as the size and quality of metabolic models increase, methods like the one presented here will be even more useful.
additional file
additional file  <dig>
figure s <dig>  comparison of runtimes for different swarm sizes.
abbreviations
apmautomatic partitioning method
biobiomass
cmcsconstrained minimal cut set
cpucentral processing unit
efmelementary flux mode
fbaflux balance analysis
fvaflux variability analysis
gagenetic algorithm
gamcsmcs software package based on a ga irrev: irreversible
lplinear program
mcsminimal cut set
milpmixed integer linear programming
pproduct
psoparticle swarm optimization
psomcsmcs software package based on pso
ramrandom-access memory
revreversible
ssubstrate
pproduct

 metabolic engineering aims to improve product yields in cellular systems by applying a variety of tools. constraint based methods which use only the stoichiometry of metabolic reactions have been particularly successful in the development of strategies towards fulfilling this aim. one important application is the prediction of knockouts to enforce desired metabolic behaviors in an organism. a method that allows one to predict efficient intervention strategies using the concept of minimal cut sets mcss, was developed by klamt and gilles. this was generalized to constrained minimal cut sets cmcs, where in addition to blocking undesired fluxes, survival of some desired fluxes is possible. the automatic partitioning method apm uses an objective function to specify the design objectives and the partitioning of fluxes into desired/undesired is done automatically to find successively larger cmcs till a global optimum is reached. previously we showed that a genetic algorithm could reach the global optimum faster than than apm. however, all these methods are applicable only to small and medium-scale metabolic networks. in a recent work by ballerstein et al., it was shown that cmcs can be directly calculated from the stoichiometric matrix. using this method, it is possible to calculate intervention strategies even in genome-scale metabolic networks. another work extended this concept to include regulation. a limitation of this method is that the desired flux or flux ratio of a metabolite has to be manually specified to get corresponding cmcs. there exist other constraint based methods for predicting intervention strategies. optknock solves a bi-level optimization problem, to predict knockouts leading to maximal product formation at maximal growth. a three-level optimization problem is used to maximize minimal product formation in robustknock. optgene uses a genetic algorithm to predict knockouts. similarly, evolutionary algorithms and simulated annealing have been used in. another metaheuristic approach was using a hybrid of bees algorithm with flux balance analysis fba. while these methods optimize for design goals, doing so with a minimal number of knockouts is not necessarily guaranteed. from an engineering perspective, we would like the organism to have a guaranteed high yield for the product of interest. given that even in the face of genetic perturbations microorganisms redirect metabolic flux towards maximizing cellular growth, this high yield must be maintained at high growth rates. additionally, the number of knockouts should be as small as possible to facilitate easy implementation in the laboratory. here we present a new method, psomcs, which uses particle swarm optimization pso along with the method developed in to calculate cmcs while overcoming the mentioned limitations of other methods. our basic motivation is to combine the computational rigour of cmcs with the flexibility of the optimization-based approaches in order to solve intervention problems efficiently. we aim to find not only the optimal intervention strategy for a given design but also the best possible design. in addition, we show that psomcs is also faster than other methods which try to find cmcs leading to optimal design objectives. methods calculating cmcs a metabolic network of m internal metabolites connected by n reactions in steady state is represented by the set of linear equations nr= where n is a m×n matrix consisting of stoichiometric coefficients of all participating reactions such that each column represents one reaction. r is a vector of reaction fluxes. reactions can be both reversible and irreversible, thereby imposing the constraint ri≥0∀i∈irrev. and define a flux space. depending on the desired outcome, an intervention problem can be set up dividing this space into desired and undesired fluxes. the set of undesired fluxes for t reactions can be defined by tr≤t where t∈ℝt×n and t∈ℝt× likewise, the set of desired fluxes for d reactions can be defined by dr≤d with d∈ℝd×n and d∈ℝd× in, cmcs are calculated by first solving a series of mixed integer linear programming milp problems representing and and then filtering those solutions which also satisfy. in, this is combined into a single system represented as in ) nrevtirev−irevtrevt0nirrtiirr−iirrtirrt00000n0000d×uvpvnwr=≥=≤000dttw≤−cu∈ℝm,vp,vn∈ℝn,d∈ℝd,vp,vn,w,rirr≥ c> note that the n and t matrices have been split into reversible and irreversible submatrices. similarly, identity submatrices for reversible and irreversible reactions are represented by the matrices i rev and i irr respectively. cmcs are directly calculated by finding solutions with minimum number of non-zero entries in v p,v n. additionally binary indicator variables z p and z n are introduced such that z p i= if v p i= and z p i= if v p i> and similarly for zn, vn. only one direction of v can be active, hence zpi+zni≤ we set up the following optimization problem minimize∑i=1ns.t., with the additional constraint that the flux through a reaction is turned off if it is part of a cmcs, i.e., r i= if z p i= || z n i= with this system it is possible to find cmcs which will result in designs satisfying constraints on yields/fluxes specified by,. however, we would like to have a method which given some design objectives calculates cmcs corresponding to optimal values for the design objectives. since any design can be represented as a function of t,d,t and d, the optimization problem can be stated as maxfs.t.. in other words, the problem is to find optimal combinations of {target/desired} yields for all reactions to be optimized. this is not easy for a few reasons. in general, this is a non-linear optimization problem. non-linear optimization is known to be inherently complex with general deterministic solutions being impossible to find. secondly, slight adjustments in, could result in completely different cmcs with different cardinalities. finally, not all such combinations will result in cmcs. these issues become acute when the search space is more dense with many possible combinations, as in large and genome-scale metabolic networks. we attack this problem using pso as it has been successfully used to find solutions to complex non-linear optimization problems in other fields. particle swarm optimization pso is a metaheuristic inspired by the flocking behavior of birds. in pso, particles distributed within a multidimensional space collectively move towards an optimum guided by a fitness function. particle fitness is determined by its position in the search space. the motion of a particle is influenced by its neighbours and the currently known fittest particle. more information on pso can be found in. typically, a particle is made up of three j-dimensional vectors, where j is the dimensionality of the search space. these represent the current position x, its previous best position p which is the position corresponding to the highest fitness achieved by the particle and the velocity v, fig. particle motion is guided by the following equations, vi=χ{vi+φ1β1+φ2β2} fig. schematic of the pso particle. a particle stores three types of information: the current values, values corresponding to its own previous best fitness and velocities corresponding to each objective xi=xi+vi i∈{.j}. g is the position corresponding to the global best fitness of the entire swarm till the current t. φ are called “acceleration constants” and determine the relative influence of the particle’s own knowledge and that of the group, both of which are commonly set to. β are uniformly generated random numbers within the range needs constraints specified by, to calculate corresponding cmcs. for example, consider a network which has, among other reactions, a substrate uptake reaction r s, a reaction for the product secretion r p and one for biomass r bio. an optimal design could be stated as having r p/r s≥x and also that biomass fluxes of r bio/r s≥x exist. however, we don’t know the combinations of x x resulting in optimal design. this is where a pso can be useful. after initializing x, v, the set of positions and velocities for all particles, within the range of values for {x x 2} on some constant r s, the pso iteratively finds increasingly better solutions for using and and moves towards the global optimum. the psomcs flowchart is shown in fig. fig. flowchart of psomcs. p and g are the current particle best and global best respectively. the algorithm stops when the number of iterations reaches a pre-specified maximum or if the maximum fitness remains unchanged for a pre-specified number of iterations the fitness function will depend on the nature of the desired optimum. considering that our objective is to have a design with high yields and minimal knockouts, the following fitness function was used, f=1−|cmcs|n·∏ixixi.  to clarify the working of psomcs, we first apply our method to a small toy network, optimizing for only a single reaction. next, to confirm the accuracy of our predictions, we compare our method against another method based on a genetic algorithm which we had previously developed. the model used is the medium-scale e. coli core model presented in. finally we find optimal intervention strategies for maximizing the minimal product yield in a genome-scale metabolic network. fba was used to calculate the range of yields for each objective and particles were initialised within this range. only one solution is calculated for a milp. the parameters used are shown in table implementation of psomcs was done using perl http://www.perl.org/. for the performance critical parts of the program, i.e., solving the milp and also the lp, the ibm ilog cplex optimization studio - a commercial optimization package - was used through the math::cplex perl module. also, our algorithm is designed to make use of modern cpu architectures and can be run in parallel on multiple cores. details of parameters used for the different models consider the network given in fig. we wish to find minimal knockouts which will ensure the highest possible yield for reaction r in the first iteration, cmcs corresponding to low yields are found. in the second iteration, all particles move towards higher yields. one particle, on the solution of its dual system gives the cmcs of ‘r r9’. removal of r and r from the network blocks all flux through r and r thus redirecting the network flux through r this corresponds to the highest minimal yield of for r fig. psomcs small example. running the psomcs on a toy network. this network has three input reactions, which can be assumed to be substrates and three secretion reactions, which can be assumed to be three different products. we want to maximise the yield of r that is maximize ). note that the particles operate in a single dimensional search space and x represents the yield for r after performing fba to determine the maximum and minimum yields for r given unit substrate uptake, four particles are initialised within this range. initial velocities are also assigned. cmcss are calculated after creating and solving the dual system. fitness is a function of x and the cardinality of the cmcs. g corresponds to x with the highest fitness which is particle after both the first and second iterations. after the first iteration, every particle except the first has a value for p. note that for particle a yield higher than is guaranteed. in reality, the minimal yield with the corresponding cmcs is which is also the case for particle this is the value the algorithm will return if allowed to run for a few more iterations we apply psomcs to generate designs in an e. coli core network which will ensure high yield of ethanol even in the face of high growth. this network was previously used to design a high yield ethanol producing strain in. this model has reactions and metabolites. we had previously used this model to predict optimal intervention strategies using a genetic algorithm, which we had shown to be faster than other current approaches, particularly compared to apm, which is guaranteed to find the optimal solution. here we compare our approach with gamcs in terms of speed and accuracy of results. the machine used had the following specifications − cpus, cores, intel xeon x ghz, running an ubuntu lts operating system. the time taken for a typical psomcs and gamcs run is plotted in fig. 4 a. the superiority of our method in terms of speed can be clearly observed. gamcs takes seconds to reach the maximum fitness. psomcs takes only seconds for the same. this is an over fold improvement in performance. in comparison, apm would not only require that the desired efms be assigned weights, but also the time taken by it would have been outside the boundaries of this plot. the cmcs corresponding to the optimum obtained by both gamcs and psomcs are exactly the same. figure 4 b is one of the designs corresponding to a high fitness. this design was in the solution pool of both the pso and ga methods. in this design, a minimum ethanol yield of is guaranteed even when the growth rate is. also, as can be expected, production of competing by-products: acetate, lactate and succinate is blocked. additionally, flux through the oxidative part of the pentose phosphate pathway is blocked and so is the pyruvate-malate cycling. multiple cmcs resulting in similar design characteristics were returned by our method. fig. comparing the runtimes of psomcs and gamcs. a plotting the runtimes of psomcs and the ga we had previously implemented clearly shows psomcs is orders of magnitude faster than gamcs. note that the time axis is logarithmic and that both algorithms reach the same maximum fitness. the change in knockout sizes are indecated by the numbers along the line. b both methods also produce similar designs, an example of which is shown. this design is obtained with knockouts. the plot was generated by applying the knockouts on the complete set of efms of the escherichia coli core model. r_norm is the sum of uptake rates for the five carbon substrates, glucose, galactose, mannose, arabinose and xylose under aerobic conditions to test the capabilities of our method we applied it to the genome-scale model of e. coli presented in. our aim was to find cmcs that result in an scenario of growth-coupled ethanol yield. a few strategies were used in to reduce the network size. these strategies are aimed at reducing the network size and improving computational efficiency, which takes real growth conditions into account and removing all superfluous components. first, the network was reduced to grow anaerobically on glucose as the only carbon source. the resulting network has reactions and metabolites. network compression was done by combining reactions operating at fixed ratios into reaction subsets. exchange reactions, spontaneous reactions and reactions essential for the ethanol and biomass production were excluded from participating in cmcs by setting their corresponding zp, zn variables to zero. the machine we used for this test had cpus, 396gb ram, intel xeon e5- ghz processor, running on ubuntu lts. the cmcs cardinality was limited to with particles being processed in parallel, the program was run for iterations. it took iterations to find the optimal design. one of the designs is shown in fig. along with designs obtained using optknock and robustknock on the same machine. the envelope of the strain specific phenotypic solution space was calculated with flux variability analysis fva of the iaf network while considering the respective knockouts predicted by each method. the minimally required biomass production was set at and both were limited by unit glucose uptake and a maximum knockout size of optknock took minutes to run while robustknock ran for minutes. the minimal ethanol yields were in both cases. as can be observed, psomcs offers a better design with the ethanol production being strongly coupled to biomass production and at no point falls below a yield of. fig. design for a genome-scale e. coli model. e. coli was designed for enhanced ethanol production using the genome-scale iaf model. for comparison, designs obtained using optknock and robustknock are also presented. the design using psomcs guarantees a minimal ethanol yield of, in contrast this is for both robustknock and optknock. all designs have a maximum biomass production rate greater than with the one for psomcs being comparatively lower. the maximum yield for all the designs is the given plots have been generated by using fva on the iaf model while considering the respective knockouts produced by each method. the fba solution space at maximum growth is highlighted, with crosses indicating the maximum and squares the minimum ethanol yield. all designs involve knockouts - for psomcs, for optknock and for robustknock discussion here we have presented a method, psomcs, to design strains with high minimal product yield using knockouts of minimal possible size. to do this, we employ a pso together with the direct enumeration of cmcs developed in. this method has made it possible to find cmcs in large and genome-scale networks. however, it is not designed to optimize engineering goals. that is, we we would like to find not only the optimal intervention strategy for a given design but the best possible design too. finding intervention strategies that achieve this is an important goal of metabolic engineering, especially in the production of industrially important chemicals. we deliver on this goal by using a pso built on top of the base provided by the direct enumeration of cmcs. our method thus expands the utility of this method. additionally we would like to point out that in the case of optimizing for a single reaction, solving with continuous values within the range for that reaction would suffice. however, in the presence of multiple objectives this task becomes computationally exhaustive and infeasible, thereby justifying the use of a metaheuristic approach such as the one used here. there have been other methods with a similar strategy as ours, which is the use of a metaheuristic in combination with another method like linear programming. most methods have relied on genetic algorithms, evolutionary algorithms and simulated annealing and also an artificial bees algorithm. ours is the first attempt at using the dual method in a similar fashion, along with the use of a pso. as shown by the comparison with optknock and robustknock in fig. although all designs have the same highest ethanol yield of psomcs provides a design with the highest guaranteed minimal ethanol yield. robustknock was developed to overcome the ’too-optimistic’ nature of optknock and this is reflected in the nature of their respective designs. also of note is the fact that both optknock and robustknock need a minimal level of biomass production to be manually specified while psomcs does not. in fact, if we reduce the minimal biomass production requirement to, robustknock runs for over hours without finding the optimum. running optknock and robustknock multiple times with different biomass levels will result in different solutions, some of which will be better than others. psomcs eliminates this need to manually set reaction fluxes and searches the entire feasible space of biomass yields to find the optimal one. growth-coupling is a key principle in metabolic engineering. it requires that growth should only be feasible if a desired compound, like ethanol, is mandatorily produced as by-product. it can be seen in fig. that psomcs achieves this with a growth rate about one third of the wild-type. however, growth-coupling does not enforce nor require that the maximal product yield is attained at a non-zero growth rate. in fact fig. illustrates the rule rather than the exception, as typically the maximum product yield is achieved at zero growth. furthermore, an ideal production state will be characterized by zero growth, where all available resources are used for product formation. in this senses, biomass production can be seen as an “unwanted” by-product. recent advances in fermentation processes employ zero-growth approaches. however, these approaches are associated with many challenges which go far beyond the scope of the presented work. nevertheless, fig. indicates that the presented designs retain their wild-type behavior to be operated as optimal zero-growth factories. in heuristic search algorithms, performance comes at the cost of being too specific to the problem being solved. by virtue of having few parameters, psos are less affected by this problem. in our implementation, we have used parameter values as found in the general pso literature without the need to adjust them. the only parameters that we adjusted were the number of particles and the number of iterations. we clearly use fewer particles than is typical. this is because we found a population size of to be sufficient for our needs. although we have sampled the entire solution space, particles can easily be forced to explore a subspace. certain reactions can be excluded from being considered for knockouts by forcing their corresponding indicator variables in the dual system to be our fitness function is specific to our target design, however new fitness functions can be thought of depending on the desired final objective. our method produces cmcs leading to designs with similar characteristics as the one used in. our method also returns multiple solutions. the limiting factor in our method is the milp for the dual system. milps are more difficult to solve than lps and may consume large amounts of time as well as memory. during our runs, the search tree generated by cplex’s branch and cut algorithm for a single milp grew to consume over gb of memory when limited to a knockout size of this memory consumption grows quickly with increasing knockout size, thereby limiting the ability of psomcs to find the optimal solution. improvements in run time can be made by forcing psomcs to explore only a part of the flux space leading to a smaller solution space to be explored. for instance lets consider the design in fig. with a minimal biomass yield of, the optimal design presented here was found within hours. further improvements to performance could be obtained by following the strategies outlined in. also, algorithmic improvements in solving milps could be useful in this regard. here we have dealt only with knockout strategies to design better strains. it can easily be extended to include the concept of regulatory mcs introduced in which combine reaction up/downregulation with knockouts. there are other constraint based methods dealing with intervention strategies like gene knock-ins and up/downregulation. psos and swarm intelligence algorithms in general may be used to compliment these methods.  psomcs finds the best possible design in metabolic networks given multiple objectives with the corresponding cmcs. we have demonstrated its capability in finding optimal knockouts and designs in genome-scale metabolic networks. it finds competitive designs compared to standard tools and is orders of magnitude faster than efm based tools in finding the optimal solution. psomcs could be used to predict minimal knockouts resulting in optimal yields in industrially important microorganisms. as the size and quality of metabolic models increase, methods like the one presented here will be even more useful. additional file additional file figure s comparison of runtimes for different swarm sizes. abbreviations apmautomatic partitioning method biobiomass cmcsconstrained minimal cut set cpucentral processing unit efmelementary flux mode fbaflux balance analysis fvaflux variability analysis gagenetic algorithm gamcsmcs software package based on a ga irrev: irreversible lplinear program mcsminimal cut set milpmixed integer linear programming pproduct psoparticle swarm optimization psomcsmcs software package based on pso ramrandom-access memory revreversible ssubstrate pproduct 
knockout strategies, particularly the concept of constrained minimal cut sets, are an important part of the arsenal of tools used in manipulating metabolic networks. psomcs finds competitive knockout strategies and designs compared to other current methods and is in some cases significantly faster. to illustrate the working of psomcs, we apply it to a toy network. psomcs not only finds solutions comparable to previously published results but also it is orders of magnitude faster. we would however like to find not only the optimal intervention strategy for a given design but the best possible design too. next we show its superiority by comparing its performance against other comparable methods on a medium sized e. coli core metabolic network. finally, we use psomcs to predict knockouts satisfying multiple objectives in a genome-scale metabolic model of e. coli and compare it with optknock and robustknock. given a specific design, cmcss can be calculated even in genome-scale networks. it can be used in identifying knockouts which will force optimal desired behaviors in large and genome scale metabolic networks. systems biologymetabolic networksdual metabolic networkminimal cut setsstrain optimizationknockoutsmetabolic pathway analysisaustrian research promotion agency ffgp issue-copyright-statement© the author 2017
 knockout strategies, particularly the concept of constrained minimal cut sets, are an important part of the arsenal of tools used in manipulating metabolic networks. given a specific design, cmcss can be calculated even in genome-scale networks. we would however like to find not only the optimal intervention strategy for a given design but the best possible design too. our solution is to use particle swarm optimization along with the direct calculation of cmcss from the stoichiometric matrix to obtain optimal designs satisfying multiple objectives.  to illustrate the working of psomcs, we apply it to a toy network. next we show its superiority by comparing its performance against other comparable methods on a medium sized e. coli core metabolic network. psomcs not only finds solutions comparable to previously published results but also it is orders of magnitude faster. finally, we use psomcs to predict knockouts satisfying multiple objectives in a genome-scale metabolic model of e. coli and compare it with optknock and robustknock.  psomcs finds competitive knockout strategies and designs compared to other current methods and is in some cases significantly faster. it can be used in identifying knockouts which will force optimal desired behaviors in large and genome scale metabolic networks. it will be even more useful as larger metabolic models of industrially relevant organisms become available. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords systems biologymetabolic networksdual metabolic networkminimal cut setsstrain optimizationknockoutsmetabolic pathway analysisaustrian research promotion agency ffgp issue-copyright-statement© the author 2017 