Optimal In Silico Target Gene Deletion through Nonlinear Programming for Genetic Engineering

Background Optimal selection of multiple regulatory genes, known as targets, for deletion to enhance or suppress the activities of downstream genes or metabolites is an important problem in genetic engineering. Such problems become more feasible to address in silico due to the availability of more realistic dynamical system models of gene regulatory and metabolic networks. The goal of the computational problem is to search for a subset of genes to knock out so that the activity of a downstream gene or a metabolite is optimized. Methodology/Principal Findings Based on discrete dynamical system modeling of gene regulatory networks, an integer programming problem is formulated for the optimal in silico target gene deletion problem. In the first result, the integer programming problem is proved to be NP-hard and equivalent to a nonlinear programming problem. In the second result, a heuristic algorithm, called GKONP, is designed to approximate the optimal solution, involving an approach to prune insignificant terms in the objective function, and the parallel differential evolution algorithm. In the third result, the effectiveness of the GKONP algorithm is demonstrated by applying it to a discrete dynamical system model of the yeast pheromone pathways. The empirical accuracy and time efficiency are assessed in comparison to an optimal, but exhaustive search strategy. Significance Although the in silico target gene deletion problem has enormous potential applications in genetic engineering, one must overcome the computational challenge due to its NP-hardness. The presented solution, which has been demonstrated to approximate the optimal solution in a practical amount of time, is among the few that address the computational challenge. In the experiment on the yeast pheromone pathways, the identified best subset of genes for deletion showed advantage over genes that were selected empirically. Once validated in vivo, the optimal target genes are expected to achieve higher genetic engineering effectiveness than a trial-and-error procedure.


Introduction
Selecting in silico, in a dynamic model of gene regulatory and metabolic networks, the right target genes for deletion so as to modify phenotypes can substantially expedite and lower the cost of genetic engineering.The target genes for deletion typically play key regulatory roles in the expression of downstream genes or metabolites to alter a phenotype to desirable states.The applications of genetic engineering are enormous.By genetically engineering plants to contain high levels of cellulose and hemicellulose [1], one may absorb the prohibitive cost of cellulose pretreatment before biomass-to-biofuel conversion.The brain tumor therapy using genetically engineered brain cells has eradicated tumors completely and affects tumor regression [2].Current in vivo genetic engineering is often by trial-and-error, and unavoidably slow and sub-optimal.The few extant in silico genetic engineering strategies are seriously hampered by the scarcity of realistic dynamic models of gene regulatory and metabolic networks.However, we anticipate a closing gap between in vivo and in silico genetic engineering as realistic computational models of networks are made increasingly available by powerful datadriven network reconstruction software from high-throughput systems biology experiments.
Recent work by Deutscher et al. [3] and Nakae et al. [4] provides multiple gene knockout solutions to optimize the concentrations of designated metabolites in static models of metabolic networks.Our work extends to dynamic models, searching the target genes in silico from any subset of genes in a gene regulatory network (GRN) for deletion to maximize the concentration of a downstream gene.Using the probabilistic Boolean network model, Faryabi et al. [5] pose an integer programming problem to maximize the benefit of a cancer patient from the treatment which intervenes the activity of a gene over time.The problem is solved using dynamic programming in optimizing a downstream gene by turning on or off only a single target gene.Based on flux balance analysis, Alper et al. [6] and Jin et al. [7] formulate a linear programming problem, to modify the metabolic pathways in wild type E. coli.They introduce the method of minimization of metabolic adjustment to revise the objective function to be quadratic for mutants.Both an exhaustive search and a greedy algorithm have been employed to optimize the yield of lycopene synthesis in the metabolic network by overexpressing or deleting three genes.They show that deleting three genes improves the phenotype of interest more effectively than deleting a single gene.
Motivated by the three-gene-deletion advantage, we study the more general multiple gene knockout (GKO) problem.Although we call all variables gene in our terminology, a variable can represent the concentration of a protein, an mRNA, or a metabolite.We use the discrete dynamical system (DDS) model to represent GRNs [8][9][10][11][12][13].DDS models can be reconstructed from observed trajectories through data-driven methods [14][15][16], some of which can run on parallel supercomputers such as [13].A nonlinear integer programming problem is formulated to define the GKO problem.We prove the nonlinear integer programming problem to be NP-hard.To approach efficiently the global maximum of the nonlinear integer programming problem with a generally non-concave objective function, we transform it to a nonlinear programming problem with fewer decision variables.We offer an algorithm called GKONP to solve the nonlinear programming problem.GKONP prunes insignificant terms in the objective function and takes advantage of the differential evolution algorithm, a parallel global optimization method.We use both the yeast pheromone pathway model and simulated models to demonstrate the performance of the GKONP algorithm.

Mathematical Formulation of the GKO Problem
We introduce the DDS model and formulate a nonlinear integer programming problem to search the optimal regulatory target genes for deletion.Here, we give the problem definition and notations.
The DDS Model.We use the DDS model [13] to represent dynamical interactions in GRNs.DDS modeling is data-driven and has been used for characterizing the cell cycle network [9].The model assumes that the change rate of each gene at the current time point is a linear combination of concentrations of genes at the previous time point.Thus state transitions are independent of each other.Let N be the number of genes.Let t be the discrete time starting from 0. Let h denote the actual time between two consecutive discrete time points.Let g i ½t be the concentration of gene i at time t.Let g½t~(g 1 ½t, g 2 ½t, . . ., g N ½t) T be a state vector of concentrations of all genes at time t.Then, the 1st-order linear DDS model is defined by where Q is an N|N regulation matrix, epitomizing a GRN.Q can be estimated with experimental data from wild type under normal and perturbed conditions.Letting A~hQzI, we have We call A the system matrix.Evidently the solution to the DDS model is Let a ij be the entry at row i and column j of matrix A. a ij is zero if gene j is not a parent (regulator) of gene i.Matrix A is sparse when the number of parents of each gene is small.: Thus, the knockout DDS solution is g½t~(diag (x)A diag (x)) t g½0.
Using the DDS solution, we define the GKO problem to maximize the objective function f (x), denoting the concentration of gene z at time T, by knocking out a subset of genes: Let x Ã be an optimal solution to the GKO problem.As we want to maximize the concentration of downstream gene z, it should not be considered for deletion and hence the constraint x z ~1.Notations -Path, Weight, and Contribution.We define path, weight of a path, and contribution of a path, to be used in the rest of the paper.A path from gene i at time 0 to gene z at time T over T time steps is a Tz1 dimensional vector, (k 0 ,k 1 , . . .,k T ) T , where k 0 ~i and k T ~z.The path is illustrated in Fig. 1.
The weight of a path is W (k 0 ,k 1 , . . .,k T )~P T t~1 a kt,kt{1 .The contribution of gene i to g z ½T through a path is defined by g i ½0 : W (k 0 ,k 1 , . . .,k T ).A path is negative/zero/positive if the contribution through the path to g z ½T is negative/zero/positive, indicating whether gene i influences g z ½T negatively or positively.

Time Complexity of the GKO Problem
We show that the GKO problem is NP-hard by reducing the NP-complete vertex cover problem to a special case of the GKO problem.Let G~(V ,E) be an undirected graph with a set V of n vertices and a set E of edges.A vertex cover is a subset of V that contains at least one end point of each edge in E. The vertex cover problem is to find a smallest vertex cover of G.We use C to represent the indices of vertices in a vertex cover of G.The Vertex Cover Problem Is a Special Case of the GKO Problem.We construct a (2nz1)|(2nz1) matrix, A, from graph G by with where A 1 is an n|n symmetric matrix, the row and column of whose non-zero entry corresponds to an edge in G, B 1 is a diagonal matrix 2I n|n , {1 v is a 1|n matrix whose entries are all {1, 1 v is a 1|n matrix whose entries are all 1, and 0 n|n The 2nz1 genes in the GKO' problem can be separated into three groups by their indices: f1, . . .,ng, fnz1, . . .,2ng, and f2nz1g.Only paths originating from group one, shown in Fig. 2, influence g 2nz1 ½2.All other paths to g 2nz1 ½2, not shown, originating from either group two or three, contribute zero to g 2nz1 ½2.We further define three types of paths, shown in Fig. 2, all originating from some genes in group one, as follows: N Type 1 path: it goes from a gene in group one at time 0, via another gene in group one at time 1, to gene 2nz1 at time 2 with a weight of 3 : {1~{3.Both gene i and j contribute {3 through their corresponding type 1 paths if (v i ,v j ) is an edge in graph G and both gene i and j exist in the network.Therefore, the number of type 1 paths is the number of nonzero elements in A 1 .
N Type 2 path: it goes from group one, via group two, to gene 2nz1 with a weight of 2 : 1~2.A gene in group one contributes 2 to g 2nz1 ½2 through its corresponding type 2 path if it exists in the network.Therefore, the number of type 2 paths is the number of existing genes in group one.
N Type 3 path: it goes from group one, via gene 2nz1, to gene 2nz1 with a weight of {1 : 1~{1.A gene in group one contributes {1 to g 2nz1 ½2 through its corresponding type 3 path if it exists in the network.Therefore, the number of type 3 paths is the number of existing genes in group one.
As the initial state is non-negative, no genes in group two should be knocked out, because doing so would not possibly increase g 2nz1 ½2, due to type 2 paths being non-negative.Thus, we consider deleting genes from group one as only feasible solutions to the GKO' problem.
Evidently, it takes polynomial time O(n 2 ) to construct the GKO' problem from the vertex cover problem.
The GKO' Problem and the Vertex Cover Problem Are Equivalent.The two problems are equivalent if and only if any smallest vertex cover C Ã of the vertex cover problem translates to an optimal solution x Ã to the GKO' problem and vice versa.
Let X be the set of all feasible solutions to the GKO' problem.Let H be the power set of f1,2, . . .,ng representing all subsets of vertices in G.We define a bijective function, w, from H to X by Function w translates any subset h[H of vertices in G to a feasible solution x[X to the GKO' problem with a corresponding objective function value g h 2nz1 ½2.When h~C Ã , the objective function value is Proof.By Fig. 2, there are three types of paths influencing g C 2nz1 ½2.Since C is a vertex cover for graph G, A 1 in equation ( 7) of matrix A(w(C)) is a zero matrix.That means there is no network between any two genes in group one and, then, genes contribute nothing to g 2nz1 ½2 through a type 1 path if we delete all gene i for all i in C from the GKO' problem.However, each gene i in group one, which is not deleted, contributes two through a type 2 path and negative one through a type 3 path.Therefore, Proof.Since C Ã is the index set for a minimum cover, we have According to Lemma 1, we have.
Let h be a non-vertex-cover subset of vertices.Let C be a smallest vertex cover that subsumes h.Then g h 2nz1 ½2ƒg C 2nz1 ½2 holds true.Proof.According to Fig. 2, one additional type 1 path contributes {3 to gene 2nz1 at time two while one additional type 2 path contributes 2 to gene 2nz1 and one additional type 3 path contributes {1 to gene 2nz1.
Since C is a vertex cover and h belongs to C, genes in group one have several additional paths to contribute nonzero values to g 2nz1 at time two if we only delete gene i for all i in h instead of in C. Let the difference of sets C and h be C D .One more gene adding into the network from C D causes more than one additional nonzero element in A 1 in equation (7).We know that the number of type 1 paths is the number of nonzero elements in A 1 .Therefore, the total contribution from the additional type 1 paths is less than As the number of type 2 or 3 paths is the number of existing genes in group one, the contribution from the additional type 2 paths is and that from the additional type 3 paths is The total contribution of those additional paths is less than Since the value in equation ( 17) is negative, value g h 2nz1 ½2 is less than g C 2nz1 ½2 and this lemma is proved.Combining Lemmas 2 and 3 establishes that g h 2nz1 ½2ƒg C Ã 2nz1 ½2 for any subset h[H of vertices in G if C Ã is a smallest vertex cover.Let x Ã be an optimal solution of GKO' and g Ã 2nz1 ½2 be its maximal value.We have the following two propositions.
. By definition of g Ã 2nz1 ½2, it is also impossible to have g C Ã 2nz1 ½2w g Ã 2nz1 ½2.Therefore, we must have g Ã 2nz1 ½2~g C Ã 2nz1 ½2.Proposition 5. Let C o be w {1 (x Ã ).Then, C o is a smallest vertex cover of G.
Proof.(By contrapositive) Assume C o ~w{1 (x Ã ) is not a smallest vertex cover of G. Then one can find a smallest vertex cover C Ã of G. Thus, it must follow by either Lemma 2 or 3 that g C Ã 2nz1 ½2wg Ã 2nz1 ½2, which contradicts the fact that g Ã 2nz1 ½2 is maximal.Therefore, C o must be a smallest vertex cover with DC o D~DC Ã D.
Propositions 4 and 5 establish that the GKO' and the vertex cover problems are equivalent.
The GKO Problem is NP-Hard.Theorem 6.The GKO problem is NP-hard.Proof.By Propositions 4 and 5, any solution to the vertex cover problem translates to a solution to the GKO' problem and vice versa.Since the vertex cover problem is in its most general form, any instance of the vertex cover problem is thus reducible to the GKO' problem.As the vertex cover problem is NP-complete and it can be reduced in polynomial time to the GKO' problem, a special case of the GKO problem, the GKO problem is NP-hard.

The Approximation Algorithm of GKONP
As the number of feasible solutions to the GKO problem increases exponentially with network size N, it is impractical to solve it by exhaustive search when N is large.Using the concept of paths, the GKO problem is rewritten to an equivalent nonlinear programming problem.Combining a strategy on pruning the insignificant terms in the objective function and a differential evolution algorithm, we provide a heuristic algorithm to the NPhard GKO problem.
Nonlinear Programming for the GKO Problem.We rewrite the nonlinear integer programming problem to an equivalent nonlinear programming problem.Let P T (i,z) be the collection of paths from gene i to gene z over T time steps.The sum of contributions of various paths from gene i to z over T time steps is It follows that the objective function f (x) of the GKO problem is the sum of contributions from all genes to gene z: A path (k 0 ,k 1 , . . .,k T ) may visit a gene more than once.We extract the unique genes on the path to form a set fk Only a negative path, (k 0 ,k 1 , . . .,k T ), in P T (i,z) gives a negative term, W (k 0 ,k 1 , . . .,k T )g i ½0, in equation (20).Therefore, we shall delete genes in negative paths to maximize equation ( 20) and those genes only on non-negative paths need not to be considered for deletion.We denote the collection of genes on negative paths to gene z by S { T (z).Then, the size of feasible solutions of the nonlinear integer programming problem is scaled down from 2 N to 2 DS { T (z)D .
Let S { T (i,z) be the collection of genes on negative paths from gene i to gene z.Let fk Proof.
, is an optimal solution but not a vertex.Therefore, there must exist an element 0vx r v1 in the solution.Then, the value of objective function (equation 21) with this solution is If the value of function (equation 22) is positive at point ), we can increase x r to one to improve the value.Otherwise, we decrease x r to zero.Since we can improve the value of objective function (equation 21) by moving x Ã to a vertex of the hypercube search space, this lemma is proved.By Lemma 7, the original GKO problem becomes Nonlinear Programming for the GKO Problem The Filter-Dynamical-Path Algorithm.We introduce the Filter-Dynamical-Path (FDP) algorithm to approximate the objective function f (x) in the form of equation ( 23).The FDP algorithm, generating the terms of objective function f (x) step by step backward from time T to time 0, discards insignificant terms at each step.Since the long run behavior of most GRNs shall be stable, A t in the DDS model also has to be such when t increases.The contributions of most paths will thus vanish and the corresponding terms are removed by the FDP algorithm when time t is long enough.
Let P z (j,t,z,T) (P { (j,t,z,T)) denote the collection of those positive (negative) paths through gene j at time t to gene z at time T and their weights.P(j,t,z,T) represents the union of P z (j,t,z,T) and P { (j,t,z,T).W z j ½t (W { j ½t) denotes the total weight of positive (negative) paths in P z (j,t,z,T) (P { (j,t,z,T)).The FDP algorithm, moving backward over time, removes those positive (negative) terms such that the total weight of the related paths is at most s of the total weight of the remaining positive (negative) paths.We call s the prune coefficient.Let f f be the approximate value to the true objective function value f .In our simulation study, the relative error is roughly bounded by The inequalities suggest that the smaller s is, the closer the approximation is to the true value.For instance, when s is 0:001 and T is 10, we have 0:99f ƒ f f ƒ1:01f .Details of the FDP algorithm is shown as Fig. 3.
The GKONP Algorithm.Based on the nonlinear formulation, we develop a heuristic algorithm to solve the GKO problem.We call it the GKONP algorithm, shown as Fig. 4.
It combines the FDP algorithm and a differential evolution (DE) algorithm for nonlinear programming.The GKONP algorithm simplifies the objective function first by the FDP algorithm and then use the DE algorithm to obtain a final solution to the GKO problem.
The DE algorithm [17,18] approaches a global maximum of non-concave objective functions as in the GKO problem.The DE algorithm is an evolutionary optimization method.The first step is to generate an initial population of feasible solutions, typically 2 to 50 times of the decision variables.Each individual in the population either remains unchanged or mutates to a new feasible solution in one iteration of evolution.The occurrence of a mutation depends on a trial vector and a probability p.The trial vector combines three other individuals, randomly chosen from the population.If the trial vector is a feasible solution and improves the value of objective function, then the individual mutates to the trial vector with probability p.Since the evolution of an individual is independent of others, evolutions of individuals can progress simultaneously and hence can be done in parallel.

Results
The GKONP algorithm is applied to improve the concentrations of downstream proteins or protein complexes Fus3PP, Fur1PP-Cdc28 (complex N) and Fur1PP-Gbc (complex M), involved in the yeast pheromone pathways.Moreover, we evaluate our algorithm on randomly generated DDS models to illustrate its empirical accuracy and running time.

Optimal Deletion in the Yeast Pheromone Pathways
We demonstrate our GKONP algorithm using a realistic Saccharomyces cerevisiae pheromone pathway model developed by Kofahl et al. [14], shown in Fig. 5.The model is obtained after they studied cell cycle arrest, mating activity, and pheromone sensitivity.The model is publicly available from the BioModels database [19] in the form of a dynamical system model composed of ordinary differential equations (ODEs).The pheromone signaling pathway involves a series of biochemical reactions starting with the receptor of MAT a receiving pheromone a factor from haploid MAT alpha .From the cytoplasm, the pheromone signal enters the nucleus to express downstream protein Fus3PP, protein complex N and protein complex M, which together control pheromone sensitivity, cell polarity and cell cycle arrest for preparation of cell fusion between two mating haploid yeast cells, MAT alpha and MAT a .Haploid MAT a cannot stop cell cycle to mate with MAT alpha if the concentrations of the three protein products are low.Therefore, it is desirable to engineer the yeast improve these downstream protein products to increase mating activity.
Thus, we applied the GKONP algorithm to identify upstream knockout genes to improve the concentrations of the three downstream protein products.By simulation using the ODE model, we first generated continuous-time trajectories.Second, we sampled them every 0.6 seconds from 0 to 6 seconds to obtain discrete-time trajectories.Then we reconstructed a DDS model (Appendix S1) from the discrete-time trajectories using a data-driven method [13] that balances goodness-of-fit and model complexity.The DDS model captures the transient dynamics in the pathway in which the three protein products are actively expressed.Using the DDS model as input, we ran the GKONP algorithm three times to search for three optimal target gene sets in the pathway for improving the concentrations of downstream products of Fus3PP, complex N and complex M, respectively.A feasible solution is any subset of {Ste2, Ste5, Ste11, Ste7, Ste20, Ste12, Fus3PP, Bar1, Far1PP, Cdc28}.The optimal target gene sets for improving each of Fus3PP, complex N and complex M are {Ste5, Ste7, Ste12},{Ste12} and {Ste12}, respectively.The optimal target genes obtained through GKONP algorithm were validated in the original ODE model.By assigning zero values to the deleted genes, we simulated the modified dynamics of the engineered ODE model.Figure 6 presents the transient dynamics, computed from the original ODE model as a validation, of the concentrations of Fus3PP, complex N and complex M in the wild type and five mutants from 0 to 6 seconds.The modified dynamics are compared with those of wild type and three observed mutants which have high concentrations of at least one of those three downstream protein products.These three observed mutants include a mutant whose G ba is overexpressed (double amount of G ba ) [20], a mutant whose Ste2 loses function (the hydrolysis of G a GTP to G a GDP is almost stopped) [21] and a mutant that has no phosphatase activity on Fus3PP (the concentration of Fus3PP is strongly increased for a long time) [22].From this figure, we demonstrate that, by deleting optimal subsets of target genes, the concentrations of all three desirable downstream protein products are higher than the wild type and the trial-and-error in vivo mutants.

Accuracy and Time Efficiency
Nine randomly generated DDS models were used to test the performance of the GKONP algorithm.The model sizes are Figure 6.The optimal solution from the GKONP algorithm vs. the wild type and three observed mutants.The GKONP optimization (top diamond lines) returned higher concentrations of desired downstream proteins (Left: Fus3PP, Middle: Complex M, and Right: Complex N) than both the wild type and the in vivo trial-and-error mutants [14].Left: the optimal target genes are Ste5, Ste7 and Ste12.Middle and right: the optimal target gene are both Ste12.The dynamics of wild type, the mutant (Sst2 lost-of-function), and the mutant without phosphatase activity on Fus3PP overlap in all three plots.doi:10.1371/journal.pone.0009331.g00610, 20 and 30, with 3 instances for each size.As each gene is only regulated directly by a very small portion of its network, without loss of generality, we constrained each DDS model such that each gene has less than three parent genes.Thus, no more than three entries were non-zero in each row of matrix A, whose values were between {1 to 1.Each gene had one unit of concentration at time 0. We applied the GKONP algorithm on the DDS models to maximize the concentration of the first gene at time 10.Prune coefficient s was set to 10 {5 .We also ran the brute force exhaustive search algorithm on the DDS models as a reference for performance evaluation of the GKONP.
Table 1 shows the GKONP algorithm approaches optimal solutions accurately because all six approximations for the 10 gene and 20 gene DDS models are identical to the optimal values, given by the exhaustive algorithm.
The time as a function of model sizes is shown in Fig. 7.The exhaustive search algorithm for the 30-gene DDS models took more than five days and we rounded the time to five days.The speedup of the GKONP algorithm ranges from 0.09 to 1.42 in the 10-gene models, 47 to 11,388 in the 20-gene models, and 4,763 to 165,390 in the 30-gene models.Therefore, Figure 7 suggests that the GKONP algorithm is much more efficient than the exhaustive search algorithm.
In Fig. 8 the number of paths decreases significantly after the FDP algorithm is applied, so that the DE algorithm runs in a reduced search space.
Since the GKONP searches negative paths of a DDS model, it was slower with 10 genes than the exhaustive search.However, as the number of genes increases to 20 and 30, our approach has extraordinary speedup over the exhaustive search.With the same model size, the running time advantage of GKONP becomes evident when the topology of a DDS model contains either few negative paths or very few genes in negative paths.For instance, a 20-gene DDS model had 8 genes in negative paths and the GKONP yielded a speedup of 11,388 versus another 20-gene model with a speedup of only 47.
This simulation study demonstrates empirically that the GKONP algorithm has achieved good accuracy in a practical amount of running time.

Discussion
We have established that the optimal in silico target gene deletion problem is challenging, by showing that a nonlinear integer programming formulation of the GKO problem based on the DDS model is NP-hard.A nonlinear programming solution is provided that combines heuristics based on the sparsity of typical GRNs and a parallel differential evolution algorithm for nonlinear programming.Multiple simultaneous gene deletion is handled in our approach, while all existing strategies delete one gene at a time.Our algorithm GKONP has shown its substantially reduced running time and comparable accuracy with the optimal solutions using exhaustive search algorithms.Demonstration of our solution on a realistic model of yeast pheromone pathways has suggested potential impact of our work.Hopefully, ideas presented in this paper will bring out potentially harder but biologically more viable computational problems for richer formulation of the target gene deletion problem, based on more complex dynamical system models of gene regulatory and metabolic networks with additional constraints on side effects.

1 , 0 n|1 2 , 0 n|n 3 ,
and 0 n|1 4 are all zero matrices.Now, we formulate the GKO' problem of 2nz1 genes, a special case of the GKO problem, as

Figure 2 .
Figure 2. Schematic diagram for the three types of paths influencing g 2nz1 ½2.doi:10.1371/journal.pone.0009331.g002 mƒT: As each element in (x k 0 ,x k 1 , . . .,x k T ) is either zero or one, we have P T j~0

Lemma 7 .
k 1 ,...,k T )[P T (i,z) If a nonlinear programming problem has objective function f (x) (equation21) and all decision variables x i , i[fk 00 bounded by ½0,1, then there exists an optimal solution which is a vertex of the feasible hypercube.

Figure 8 .
Figure 8.Average path reduction by the FDP algorithm.The red line with circles represents the number of paths before FDP reduction while the blue line with diamonds represents the number after FDP reduction.doi:10.1371/journal.pone.0009331.g008

Figure 7 .
Figure 7. Average running time of the GKONP algorithm and the exhaustive search algorithm.The GKONP algorithm is represented by the blue line with diamonds while the exhaustive search algorithm is by red line with circles.doi:10.1371/journal.pone.0009331.g007 Optimal Target Gene Deletion through Nonlinear Integer Programming.Based on the DDS model, a nonlinear integer programming is formulated to maximize a downstream gene by searching regulatory target genes for deletion.We define the binary knockout vector x~(x 1 , . . .,x N ) T .x i [f0,1g is 1 if gene x i is intact; x i is 0 if gene i is deleted, equivalent to setting all entries on either row i or column i in system matrix A to zero.