Metabolite optimization method based on binary quantum particle swarm optimization algorithm and flux balance analysis hybrid algorithm

Abstract

Metabolite optimization method based on a hybrid algorithm of binary quantum particle swarm optimization algorithm and flux balance analysis. The method of the present invention combines the advantages of flux balance analysis FBA and improved binary quantum behavior particle swarm optimization EBQPSO algorithm to determine the optimal knockout strategy for the desired metabolites. The EBQPSO algorithm introduces the concept of neighborhood to increase the diversity of the particle swarm and promote the identification of the global optimal solution. We tested our algorithm on two genome-scale metabolic models of Escherichia coli iAF1260 and iJR904, and the results showed that the performance of EBQPSO is competitive compared with other algorithms. In addition, the encoding method used in our algorithm is rarely used in the optimization problem of GSMMs, which is an important contribution to research in this field. The algorithm we proposed has the potential to significantly improve the optimization of microbial cell factories and meet the growing demand for microbial manufacturing products in different industries.

Description

Background Art

[0001] The present application belongs to the field of metabolic engineering, and specifically relates to a metabolite optimization method based on a hybrid algorithm of binary quantum particle swarm optimization algorithm and flux balance analysis. Technical Field

[0002] The goal of metabolic engineering is to modify and optimize metabolic pathways to achieve desired cellular phenotypes, such as the production of specific metabolites. Achieving high yields and high production efficiency is a challenging task that requires a deep understanding of the complex metabolic pathways and gene regulation of the microorganisms used for production. In recent years, the demand for microbial production has grown rapidly, and there is an urgent need to develop efficient and effective methods to optimize metabolic pathways to achieve higher yields and production efficiency.

[0003] With the rapid development of gene sequencing and high-throughput experimental technology, it has become possible to reconstruct metabolic networks. Genome-scale metabolic models (GSMMs) are mathematical models reconstructed from annotated genome sequences to characterize gene-protein-reaction relationships. A variety of analytical methods have been developed for GSMMs, such as flux balance analysis (FBA), metabolic regulation minimization (MOMA), flux variability analysis (FVA), etc. These methods use computer simulations to predict the corresponding phenotypes produced by genetic intervention.

[0004] Although these computational methods are very effective in predicting the impact of a given genetic intervention on microbial metabolism, they are limited in their ability to identify the optimal genetic interference strategy that leads to the desired microbial strain. To address this issue, several two-level linear programming (LP) tools have been developed, such as OptKnock, OptStrain, and RobustKnock. However, these two-level LP methods have certain limitations, including suboptimal solutions, long computational times, and difficulty in scaling to larger systems. In addition, these frameworks require precise mathematical transformations that increase the complexity of the underlying computations and limit the maximum number of interventions allowed per simulation. In addition, these algorithms only support the application of linear constraints and objective functions, which may not accurately represent the complexity of biological goals.

[0005] To overcome these limitations, intelligent optimization algorithms, such as particle swarm optimization (PSO), have been developed to find optimal strain improvement strategies. PSO uses a swarm of particles to explore the search space and locate the optimal solution through iterative calculations. Variants of PSO, such as quantum-behaved particle swarm optimization (QPSO) and binary quantum-behaved particle swarm optimization (BQPSO), have also been developed to improve their search performance. These intelligent optimization algorithms have been integrated with constraint-based metabolic flux analysis methods to identify optimal microbial mutant phenotypes with complex and nonlinear engineering objectives as a fitness function. For example, the PSOMOMA algorithm combines PSO and MOMA to identify gene knockout sequences that can optimize ethanol production and growth rate in Escherichia coli. These impressive results demonstrate the potential of using intelligent optimization algorithms to improve microbial bioprocesses in a sustainable manner.

[0006] In this paper, we propose a new hybrid optimization strategy called EBQPSOFBA, which combines the enhanced binary quantum-behaved particle swarm optimization (EBQPSO) algorithm with the FBA method. Our algorithm avoids local optimality by introducing the concept of neighborhood to calculate local attractors and increase the diversity of the particle swarm. The encoding method used in EBQPSOFBA is applicable to large-scale GSMMs and significantly reduces the search space. In addition, EBQPSO can also be combined with other constraint-based methods to determine the optimal strain design. Our study shows that EBQPSOFBA can effectively identify the best gene knockout strategies for various product yields and focus on improving the yield of Escherichia coli. In addition, the proposed algorithm performs well or is comparable to existing algorithms in terms of effect. We also make an important contribution to this research field by highlighting the underutilized encoding methods used in our algorithm. Summary of the invention

[0007] The present invention aims to provide a strain transformation strategy optimization method based on the Enhanced Binary Quantum Particle Swarm Optimization Algorithm (EBQPSO), which is an improvement on the original Binary Quantum Particle Swarm Optimization Algorithm (BQPSO) and introduces the concept of neighborhood to calculate local attractors. Its main goal is to increase the yield of desired metabolites in the production process.

[0008] The technical solution of the present invention:

[0009] The metabolite optimization method based on the binary quantum particle swarm optimization algorithm and the flux balance analysis hybrid algorithm has the following steps:

[0010] The purpose of preprocessing is to eliminate selective reactions that may lead to undesirable results or increase the burden on the algorithm from the total set of reactions in the model, thereby obtaining a smaller set of effective knockout reactions. The preprocessing steps are as follows: First, basic reactions required for growth are removed from the list of potential knockout reactions. The next step is to exclude reactions that are not related to genes because they are usually not biologically meaningful. Finally, biomass, ATP maintenance (ATPM), and target reactions are not considered in the search space.

[0011] Step 1. Initialize parameters: set the population size, particle dimension, maximum number of iterations, shrink-expansion coefficient, neighborhood size, and number of runs to meet the needs of the optimization process.

[0012] Step 2: Generate the initial population: randomly generate a certain number of particles, each particle represents a solution to the strain transformation strategy.

[0013] Step 3: Evaluate fitness: The fitness of each particle is calculated using flux balance analysis (FBA) to measure the performance of the particle in improving its ability to produce metabolites.

[0014] Step 4: Determine the neighborhood relationship: Define a neighborhood for each particle according to the preset neighborhood size parameters, and use a ring topology as the neighbor topology of the particle, so that other particles in its neighborhood can be determined for each particle in each iteration.

[0015] Step 5: Particle position update: Solve the local attractor by calculating the best particle in the neighborhood.

[0016] (a) The particle with the best fitness is selected from the neighborhood of each particle as the best particle in the neighborhood.

[0017] (b) When updating the particle position, a new variable G is introduced, which represents the influence of another neighboring particle that the current particle may be affected by. The calculation method of G is as follows: randomly select a different particle as a candidate solution of G. If the individual optimal fitness value of the particle is better than the individual optimal fitness value of the current particle, the individual optimal position of the particle is assigned to G. Otherwise, the global optimal position is assigned to G.

[0018] (c) The local attractor is obtained by introducing the best particle in the neighborhood through multi-point crossover. Specifically, the value of each dimension of the local attractor is calculated through a loop. In each loop, a random number is first generated. According to the size of the random number, the corresponding dimension value of the individual best position, the new variable G, or the best particle in the neighborhood is randomly selected with equal probability as the value of the corresponding dimension of the local attractor. If the random number is less than 1 / 3, the individual best position is selected; if the random number is between 1 / 3 and 2 / 3, the new variable G is selected; if the random number is greater than or equal to 2 / 3, the best particle in the neighborhood is selected. After the loop is completed, a local attractor composed of all local attractor dimension values ​​is returned.

[0019] (d) Use the BQPSO algorithm to calculate the probability value for each dimension. Perform the following operations on each bit of the substring corresponding to each decision variable of the local attractor: If the calculated probability value is greater than the random number, then invert the bit. In other words, if the bit was originally 1, set it to 0; vice versa. And assign the substring to the corresponding substring of the updated particle.

[0020] (e) After completion, all updated particle substrings are combined to obtain the new position.

[0021] (f) Ensure that the updated particle positions satisfy specific constraints to maintain the validity of the problem solution space.

[0022] (g) Use the FBA method to update the fitness value of the particle, and update the individual optimal solution and the global optimal solution.

[0023] Step 6, iterative optimization: Repeat steps 3 to 5 until the preset maximum number of iterations is reached to find the best strain transformation strategy.

[0024] Step 7, output optimization results: output the global best particle (i.e., the particle with the highest fitness) as the optimal strain transformation strategy to guide strain transformation in the actual production process.

[0025] The beneficial effects of the present invention are as follows: the improved EBQPSO algorithm is used to optimize the strain transformation strategy, and the fitness of the particles is evaluated in combination with the FBA technology. The yield of the desired metabolite is taken as the goal, and the purpose of finding the best strain transformation strategy is achieved. In the present invention, the EBQPSO algorithm, as an enhanced binary quantum particle swarm optimization algorithm, has higher search capabilities and convergence speed in the process of finding the optimal solution compared to the traditional BQPSO algorithm. The introduction of the concept of neighborhood to calculate local attractors enables the algorithm to achieve a balance between global search and local search, effectively avoiding the dilemma of falling into the local optimal solution. Therefore, the present invention has higher accuracy and stability in finding the best strain transformation strategy. In short, the present invention provides a strain transformation strategy optimization method based on the EBQPSO algorithm, which is innovative and practical. By combining the FBA technology and the EBQPSO algorithm, with the yield of the desired metabolite as the optimization goal, it provides strong theoretical support and guidance for strain transformation in the actual production process. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 is the neighborhood ring topology of particle i.

[0027] Figure 2 To decode and evaluate the solution process. DETAILED DESCRIPTION

[0028] 1.1 Flow Balance Analysis (FBA)

[0029] Schilling et al. proposed the FBA method in 1999 [Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis? . Nat Biotechnol 28, 245–248 (2010). https: / / doi.org / 10.1038 / nbt.1614], which mathematically fits metabolic flux into the metabolic network. Assuming that the intracellular metabolic reactions are in a pseudo-steady state, FBA uses metabolic goals (such as maximizing biomass or minimizing adenosine triphosphate (ATP) production) as the objective function to predict the flux distribution in the entire metabolic network. For a metabolic network with metabolites and chemical reactions, the established metabolic balance equation is as follows.

[0030] Maximize / MinimizeC (1)

[0031]

[0032] v min ≤v y ≤v max (3)

[0033] Where S mn represents the stoichiometric coefficient of the mth metabolite in the nth reaction, v n represents the flux of the nth reaction, C represents the optimization target, and v min represents the lower limit of the reaction flux, v max represents the upper limit of the reaction flux.

[0034] 1.2 Quantum Particle Swarm Optimization (QPSO)

[0035] In QPSO, the population with particles in D-dimensional space is called population X. The position vector of particle i in generation t is represented by x i (t) = (x i1 (t),x i2 (t),…,x iD (t)), which is a potential solution to the given problem. The vector pbest i =(pbest i1 ,pbest i2 ,…,pbest iD ) records the individual best position (pbest) of particle i. Vector gbest = (gbest1, gbest2, …, gbest D ) represents the global best position (gbest) of all particles in the swarm. The equation of the QPSO algorithm is first shown below.

[0036]

[0037] p id =r1×pbest id +(1-r1)×gbest d (5)

[0038]

[0039] where i = 1, 2, ..., N; d = 1, 2, ..., D. Parameters r1 and r2 are two random numbers with uniform distribution in (0, 1), and mbest represents the mean best position, which is the individual best position pbest of all particles calculated by equation (4) i The average value of mbest d represents the dth dimension of mbest. The parameter α is called the shrinkage-expansion coefficient, and it is used to control the convergence speed of the algorithm. id is the local attractor p of particle i i The dth dimension of is derived from equation (5) between pbest i A random position between gbest and gbest. Update the position vector x of particle i in generation t+1 according to equation (6): id (t+1), |*| represents the Euclidean distance.

[0040] 1.3 Binary Quantum Particle Swarm Optimization (BQPSO)

[0041] In BQPSO, |*| represents the Hamming distance instead of the Euclidean distance because the particle positions are represented by binary strings. The Hamming distance is the number of different bits between two binary strings obtained by comparing the two strings.

[0042] |XY|=d H (XY) (7)

[0043] where X and Y are two binary strings, such as positions. H Get the Hamming distance between X and Y.

[0044] mbest j represents the jth bit of mbest, which is determined by the jth bit state of pbest of all particles. The specific calculation method is as follows: if the mbest obtained from equation (1) j If greater than 0.5, then mbest j =1, indicating that the jth bit of the mbest of more particles is 1; if mbest j <0.5, then mbest j = 0, indicating that the jth position of pbest of more particles is 0; otherwise, if mbest j If it is exactly equal to 0.5, then it will randomly take 0 or 1 with a probability of 0.5.

[0045] The local attractor p of particle i i It is obtained through crossover operation, including single-point crossover or multi-point crossover. First, gbest i Perform a crossover operation with gbest to generate two offspring. Then randomly select one of them as p with a probability of 0.5. i Consider the iterative equation (6) and transform it into:

[0046]

[0047] where d represents the dth dimension, corresponding to the binary substring of the dth decision variable. Therefore, x id (t+1) and p id Represents the position vector x i (t+1) and the local attractor p i The binary substring of the d-th decision variable.

[0048] The equation can be expressed in terms of Hamming distance as:

[0049]

[0050]

[0051] The function used The reason is that the Hamming distance must be an integer. The parameter r is a random number with uniform distribution in the range (0,1).

[0052] According to the following equation, Pr d , by adding p id Each bit of Pr d The probability of mutation calculation is used to obtain the new substring x id (t+1). id All dimensions of (t+1)(d=1,2,…,D) are combined into a new x i (t+1).

[0053]

[0054] Among them l d represents the length of the dth dimension, that is, the length of the binary substring of the dth decision variable. The sum of the lengths of the binary substrings of all decision variables is represented as l.

[0055]

[0056] Based on the above definition and the modification of the iterative equation, the BQPSO algorithm is described as follows.

[0057] Step 1. Initialize the binary bit arrays of all particles: X, pbest = X;

[0058] Step 2. Calculate the objective function value of each particle and compare it with the previous objective function value of pbest. If the current objective function value is better than the previous objective function value, update pbest;

[0059] Step 3. Determine the current global best position gbest among all individual best positions of particles;

[0060] Step 4. Calculate the average optimal position mbest and local attractor p for each particle i ;

[0061] Step 5. Calculate the mutation probability Pr of each dimension according to equation (8) d , and update the new substring. The specific operation is as follows: id Each bit of Pr performs the following operations: d If the bit is greater than the random number rand(), the bit is inverted. That is, if the bit was originally 1, it is set to 0, and vice versa. id Assign to the updated substring x id. Then all substrings are combined to get the new position x i .

[0062] Step 6. Repeat steps 2 to 5 until the stopping criterion is met or the specified maximum number of iterations is reached.

[0063] 2.1 Enhanced Binary Quantum Particle Swarm Optimization Algorithm (EBQPSO)

[0064] The enhanced binary quantum particle swarm optimization algorithm (EBQPSO) draws on the ideas of Nasir et al. [MdNasir, SwagatamDas, Dipankar Maity, et al. A dynamic neighborhood learning based particle swarm optimizer for global numerical optimization [J]. Information Sciences, 2012, 209 (1): 16-36.] and adds the concept of neighborhood when calculating local gravity. The introduction of the concept of neighborhood in the EBQPSO algorithm enables particles to explore beyond their local search space and better determine the global optimal solution. By defining a neighborhood for each particle, the algorithm can consider not only the position of the particle itself, but also the positions of its neighbors, which helps prevent the algorithm from falling into the local optimal value. Specifically, we use the ring topology as the neighbor topology of the particle and assume that the neighborhood size is 5. Considering N+1≡1, if the index of a particle is i, the particles with indices i-2, i-1, i+1 and i+2 constitute the neighborhood of particle i, such as Figure 1 The particle with the best fitness value in the particle neighborhood is called the best in the neighborhood, denoted by nbest. i The introduction of the neighborhood concept is a key improvement in our algorithm, which enables particles to better explore the search space, improves the efficiency of determining the optimal solution, and enhances the performance of the algorithm.

[0065] As mentioned before, in BQPSO, the local attractor p i is the best individual position pbest of particle i i The global best position pbest of all particles in the group is cross-operated. In the proposed EBQPSO algorithm, in order to update the position of particles, we will introduce a new variable G, which represents the influence of another neighbor particle that the current particle i may be affected by. The calculation method of G is as follows: randomly select a particle k (k≠i) as a candidate solution of G. If the individual optimal fitness value pbest of particle k is k The individual optimal fitness value pbest that is better than particle i i , then G is set to pbestk Otherwise, G is set to gbest. In addition, we also introduce the neighborhood optimal particle nbest by multi-point crossover i To obtain the local attractor p i The specific operation method is to calculate p by looping l times i The value p of each dimension d id In each cycle, the algorithm first generates a random number r. According to the size of r, the algorithm will randomly select pbest with equal probability. i , G or nbest i The corresponding dimension value of is taken as p id If r is less than 1 / 3, select pbest i ; If r is between 1 / 3 and 2 / 3, select G; if r is greater than or equal to 2 / 3, select nbest i After the loop is completed, the algorithm returns the id The value of p i Next, for each dimension, the probability value Pr is calculated according to the BQPSO algorithm. d . The local attractor p i The binary substring of the d-th decision variable is denoted as p id , for the substring p id Each bit of Pr performs the following operations: d If the bit is greater than the random number rand(), the bit is inverted. That is, if the bit was originally 1, it is set to 0, and vice versa. id Assign to the updated substring x id . Then all substrings are combined to obtain the new position xi. Finally, the fitness value of the particle is updated, and the individual optimal solution (pbest) and the global optimal solution (gbest) are updated. The cycle continues until the termination condition is met. The entire algorithm aims to find the optimal solution to the problem by continuously iterating and updating the position of the particle.

[0066] 2.2 Improved hybrid binary quantum particle swarm optimization algorithm and flux balance analysis (EBQPSOFBA)

[0067] EBQPSOFBA is a combination of EBQPSO and FBA, where FBA is used to evaluate the fitness value. A brief overview of EBQPSOFBA is as follows:

[0068] The algorithm first randomly generates an initial binary bit array for each particle. Next, the individual optimal solution (pbest), global optimal solution (gbest), average optimal solution (mbest), and neighborhood optimal solution (nbest) of each particle are initialized, and their fitness values ​​are evaluated using the FBA method. The algorithm then enters a loop until the termination condition is met. During the loop, the algorithm calculates the average optimal solution (mbest) and updates the neighborhood optimal solution (nbest) of each particle. For each particle, the algorithm calculates its neighborhood optimal solution, randomly selects a different particle to obtain G, and calculates the local attractor p according to the algorithm EBQPSO. i Next, for each dimension, the probability value Pr is calculated according to the BQPSO algorithm. d . Local attractor p i The binary substring of the d-th decision variable is denoted as p id , for the substring p id Each bit of Pr performs the following operations: d If the bit is greater than the random number rand(), the bit is inverted. That is, if the bit was originally 1, it is set to 0, and vice versa. id Assign to the updated substring x id . Then all substrings are combined to get the new position x i Finally, the FBA method is used to update the fitness value of the particle, and the individual optimal solution (pbest) and the global optimal solution (gbest) are updated. The cycle continues until the termination condition is met. The entire algorithm aims to find the optimal solution to the problem by continuously iterating and updating the position of the particles.

[0069] 2.2.1 Preprocessing

[0070] The purpose of preprocessing is to eliminate selective reactions that may lead to undesirable results or increase the burden on the algorithm from the total set of reactions in the model, thereby obtaining a smaller set of effective knockout reactions. The preprocessing step refers to the method of Sarma [Mutturi S.FOCuS: a metaheuristic algorithm for computing knockouts from genome-scale models for strain optimization. Mol Biosyst. 2017Jun 27; 13(7): 1355-1363. doi: 10.1039 / c7mb00204a. PMID: 28530276.]. First, basic reactions required for growth are removed from the list of potential knockout reactions. The next step is to exclude reactions that are not related to genes because they are usually not biologically meaningful. Finally, biomass, ATP maintenance (ATPM), and target reactions are not considered in the search space.

[0071] 2.2.2 Initialization

[0072] Most published papers encode solutions using a binary string equal to the number of genes or reactions in the metabolic model, where binary bit 1 indicates knockout and binary bit 0 indicates maintenance. However, for more complex metabolic network models, this encoding method may lead to a large solution space. In practice, we usually only allow a few genes or reactions to be knocked out, resulting in a large number of invalid solutions in the solution space.

[0073] Therefore, we adopt a binary encoding method to reduce the solution space, in which only the reaction knockout is represented in the solution. A binary substring is used to represent the reaction to be knocked out, and then the substrings are combined to represent all the reactions to be knocked out. Taking the E. coli iJR904 model (containing 1075 reactions) as an example, suppose 5 reactions are to be knocked out. Therefore, a reaction to be knocked out can be represented by an 11-bit binary substring, and a solution to knock out 5 reactions can be represented by 55 bits of binary. That is, l d = 11; D = 5. In contrast, using the previous encoding method requires a binary string of length 1075 to represent the solution.

[0074] 2.2.3 Fitness Evaluation

[0075] To evaluate the solution, the binary substrings of the particles in each dimension are first decoded into their decimal equivalents. The corresponding reaction fluxes are then forced to zero, thereby adding new constraints to the metabolic model. This process is performed by simulating mutants. In this work, the set of flux values ​​for all reactions is output using FBA (see 2.1 Flux Balance Analysis), but other methods can also be considered. According to the suggestions of Patil et al. [K. Patil, I. Rocha, J. J. Nielsen, Evolutionary programming as a platform for in silico metabolic engineering, BMC Bioinformatics. 6 (2005) 308. doi: 10.1186 / 1471-2105-6-308.] In order to ensure that the microorganisms can survive after knockout, only the growth rate greater than or equal to 0.1 hr was considered -1 The mutants are designed to overproduce metabolites. Next, fitness values ​​are calculated using the partially obtained flux values, which can be appropriate objectives or fitness functions. In this work, both product rate and biomass-product coupled yield (BPCY) are used as objective functions to measure the quality of the solution. The formula for calculating BPCY is as follows.

[0076]

[0077] Where mmol is millimole, hr is hour, and gm is gram.

[0078] The whole process of decoding and evaluating the solution is as follows Figure 2 shown.

[0079] Solution decoding and evaluation process. The solution (individual) is encoded using a binary representation, where the reaction to be removed is represented in the corresponding binary. Each individual is decoded by imposing additional constraints on the original metabolic model. FBA is a mathematical simulation of cell behavior, and its results can be used to calculate fitness values.

[0080] 2.2.4 Post-processing

[0081] The purpose of post-processing is to retain the necessary knockouts, i.e., simplify the optimal solution by removing all the deletions that keep the objective function unchanged. If some knockouts are not conducive to improving the objective function, then these knockouts are considered unnecessary, as implementing the knockouts is time-consuming and expensive.

[0082] 3 Experiments

[0083] 3.1 Dataset

[0084] In this paper, two E. coli genome-scale metabolic models (GSMMs), iAF1260 and iJR904, were used to evaluate the performance of EBQPSOFBA, and the properties of these two models are shown in Table 1. Both models are publicly available from the Biochemical, Genetic and Genomic (BiGG) database [Schellenberger, J.; Park, JO; Conrad, TM; Palsson, BTBiGG: A Biochemical Genetic and Genomic knowledgebase of large scale metabolicreconstructions. BMC Bioinform. 2010, 11.] and are widely used in the fields of biotechnology and microbiology research to improve the production of desired metabolites. The simulation experiments were performed in MATLAB R2016b environment using the Constrained-Base Reconstruction and Analysis (Cobra) toolbox [J. Schellenberger, R. Que, RMT Feming, I. Thiele, JD Orth, AM Feist, DC Zielinski, A. Bordbar, N E Lewis, S. Rahmanian, J. Kang, D R Hyduke, B. Palsson, Quantitative prediction of cellular metabolism with constraint-based models: The COBRA Toolbox v2.0, Nat. Protoc. 6 (2011) 1290-1307. doi: 10.1038 / nprot.2011.308.]. The FBA was evaluated using the academic license version of Gurobi optimizer 8.11. The input settings for all simulations were as follows (flux units: mmol / h). For iAF1260, glucose uptake was set to 10; oxygen uptake was set to 18.5; and ATP maintenance was set to 8.39. For iJR904, the inputs were: glucose uptake (10); oxygen uptake (18.5 under aerobic conditions and 0 under anaerobic conditions); and ATP maintenance (7.6).

[0085] Table 1: Properties of the genome-scale metabolic models used in the experiments

[0086]

[0087]

[21] Feist, AM, et al. (2007). A genome-scale metabolic reconstruction for Escherichia coli K-12MG1655 that accounts for 1260ORFs and thermodynamic information. Molecular systems biology, 3(1), 121.

[0088]

[22] Reed JL, Vo TD, Schilling CH, Palsson BO: An expanded genome-scalemodel of Escherichia coli K-12(iJR904 GSM / GPR). Genome Biology 2003,4:R54.

[0089] 3.2 Parameter settings

[0090] According to the experiments conducted, the optimal parameters of the EBQPSO algorithm are as follows: population size is 100, shrink-expansion coefficient is 1.4, and maximum number of iterations is 300. For fair comparison, the number of knockouts is kept the same as the FOCuS algorithm, which is 3, 5, and 10, respectively, except for the post-processing step. The experiment is repeated 20 times, and only the best knockout strategy is shown. The detailed parameter settings of the EBQPSOFBA algorithm are shown in Table 2.

[0091] Table 2: Parameter settings

[0092]

[0093] 3.3 Results and Discussion

[0094] Case Study 1: Comparison with the Bi-Level Optimization Algorithm

[0095] In the first case study, the performance of the proposed EBQPSOFBA algorithm was compared with two classic two-level optimization algorithms, OptKnock and ReacKnock. OptKnock is the first two-level LP programming framework proposed by Anthony et al. ReacKnock, proposed by Xu et al., is also a classic two-level programming algorithm, and its reported results for predicting the optimal gene knockout strategy are better than OptKnock. In order to evaluate the performance of the EBQPSO algorithm, the iAF1260 model was used to predict the optimal gene knockout strategy (the number of knockout genes is 5) for five product yields (i.e., D-lactic acid, succinic acid, acetic acid, glycolic acid, and formic acid) as the fitness value. The results obtained in this work and the comparison with other two-level linear programming algorithms are shown in Table 3.

[0096] As can be seen from Table 3, the target product yields obtained by EBQPSO were either equal to those of ReacKnock or higher in the case of acetic acid and formic acid, except for succinic acid. In all cases, EBQPSO produced higher target product yields than OptKnock (see Table 3).

[0097] Table 3: Comparison of EBQPSO results with other algorithms to maximize product flux

[0098]

[0099]

[0100]

[0101] a The numbers in brackets indicate the maximum production rates reported. (Unit: mmolh -1 wxya -1 )

[0102] b Unit:mmolh -1 wxya -1

[0103] c Unit: hr -1

[0104] 【7】Burgard, AP; Pharkya, P.; Maranas, CDOptKnock: A bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol. Bioeng. 2003, 84, 647–657.

[0105] 【9】Xu Z, Zheng P, Sun J, et al. ReacKnock: identifying reaction deletionstrategies for microbial strain optimization based on genome-scale metabolic network[J]. PloS one, 2013, 8(12):e72150.

[0106] Case Study 2: Comparison with other heuristic algorithms

[0107] In this case study, the EBQPSO algorithm was compared with two other heuristic optimization algorithms, namely, the ensemble-based evolutionary algorithm (SEA) or simulated annealing (SA) process (SEA / SA) and the floral pollination combined with clonal selection algorithm (FOCuS). The performance of the EBQPSO algorithm was evaluated with the goal of improving the biomass-product coupled yield (BPCY) of two products in the iJR904 model, namely, D-lactic acid (aerobic and anaerobic conditions) and succinate (aerobic conditions). SEA / SA was proposed by Rocha et al. with a variable number of knockouts. Similar to the FOCuS algorithm, the proposed EBQPSO requires the number of knockouts to be determined in advance. In this work, the number of reaction knockouts was kept the same as the FOCuS algorithm. The comparison results are shown in Table 4.

[0108] Table 4: Comparison results between EBQPSO and heuristic algorithm in terms of biomass-product coupled yield (BPCY)

[0109]

[0110]

[0111]

[0112] a The numbers in brackets indicate the maximum production rates reported. (Unit: mmolh -1 wxya -1 )

[0113] b Unit:mmolh -1 wxya -1

[0114] c Unit: hr -1

[0115]

[16] Mutturi S.FOCuS: a metaheuristic algorithm for computing knockouts from genome-scale models for strain optimization.Mol Biosyst.2017Jun 27;13(7):1355-1363.doi:10.1039 / c7mb00204a.PMID:28530276.

[0116]

[26] Rocha M, Maia P, Mendes R, et al. Natural computation meta-heuristics for the in silico optimization of microbial strains[J]. BMC bioinformatics, 2008, 9(1):1-16.

[0117] As can be seen from Table 4, except for D-lactate under aerobic conditions, EBQPSO obtained better or at least comparable results than FOCuS and SEA / SA in the other case studies. For D-lactate production under anaerobic conditions, the EBQPSO algorithm obtained the same maximum BPCY as SEA / SA and FOCuS, while the suggested knockout strategies were slightly different. All three algorithms suggested knockouts of ADHEr and GART, with the difference that EBQPSO suggested knockout of NADH8, while both SEA / SA and FOCuS suggested knockout of FRD3. For propionate production, EBQPSO (0.3511) was higher than FOCuS (0.32981), but slightly lower than SEA / SA (0.35785). However, it is worth noting that the number of knockouts for EBQPSO was only 5, while that for SEA / SA was 15. As can be seen from Table 4, the knockout lists suggested by EBQPSO and FOCuS have two common reactions, “FUM” and “GND”, but the other three reactions are different. Similarly, the knockout lists suggested by EBQPSO and SEA / SA also have two common reactions, “FUM” and “GHMT2r”. This further proves the effectiveness of the knockout strategy given by the EBQPSO algorithm.

[0118] In addition, the SEA / SA algorithm adds a greedy algorithm to help get better results, but at the cost of a lot of calculations. The FOCuS algorithm introduces additional external calculations to divide the search range through multiple filtering, thereby narrowing the search range to less than 400. The EBQPSO proposed in this paper cleverly reduces the search space through binary coding and does not add complex strategies, making the algorithm implementation simpler while achieving satisfactory results.

[0119] The increasing demand for microbial manufacturing products in various industries has led to the development of metabolic engineering to find the best strain design that can improve the yield of desired metabolites. However, predicting the impact of genetic modifications on phenotypes is a challenging task. Although researchers have developed several bilevel LP tools and intelligent optimization algorithms to help find the optimal strain modification strategy, these algorithms still have some limitations. This paper introduces an EBQPSO algorithm that incorporates the neighborhood concept into the original BQPSO. Subsequently, the algorithm is combined with FBA to calculate the fitness value to form a hybrid algorithm named EBQPSOFBA for identifying the optimal knockout strategy. The performance of EBQPSOFBA is evaluated through experiments on two GSMMs of Escherichia coli (iAF1260 and iJR904) and compared with existing bilevel linear programming and heuristic methods. In Case 1, the EBQPSO algorithm is compared with the bilevel optimization algorithm, and the yields obtained by the EBQPSO algorithm are higher or at least equal to those of the competitors for all five different products. In case 2, compared with other heuristic algorithms, the EBQPSO algorithm produced higher or at least equal BPCY scores than other algorithms for most products, except for d-lactic acid under aerobic conditions. In addition, the EBQPSO algorithm performed very well for propionate under aerobic conditions while requiring fewer knockouts. Our results show that the EBQPSO algorithm has a strong competitive performance compared with other heuristic and bilevel optimization algorithms. Therefore, the EBQPSO algorithm can be considered as an effective optimization method for improving bioproduction of chemicals because it exhibits higher or at least comparable product yields and biomass-product coupled yield (BPCY) scores.

[0120] One limitation of the EBQPSO algorithm is that it is limited to solving single-objective problems. However, in the context of engineering a desired strain, product yield and microbial growth are often considered important factors. Therefore, future research could focus on extending the algorithm to handle multi-objective problems. In addition, this study only predicted knockout strategies. In subsequent studies, we aim to expand the application scope of BQPSOFBA by predicting overexpression and heterologous expression targets, thereby providing a more comprehensive solution for strain engineering.

Claims

1. A metabolite optimization method based on a hybrid algorithm of binary quantum particle swarm optimization algorithm and flux balance analysis, characterized in that: Here are the steps: Step 1, initialization parameters: set the population size, particle dimension, maximum number of iterations, shrink-expansion coefficient, neighborhood size, and number of runs to meet the needs of the optimization process; Step 2: Generate the initial population: randomly generate a certain number of particles, each particle represents a solution to the strain transformation strategy; Step 3: Evaluate fitness: Calculate the fitness of each particle using metabolic flux balance analysis to measure the performance of the particle in improving metabolite production capacity; Step 4, determine the neighborhood relationship: define a neighborhood for each particle according to the preset neighborhood size parameter, and use the ring topology as the neighbor topology of the particle, so as to determine other particles in its neighborhood for each particle in each iteration; Step 5, particle position update: solve the local attractor by calculating the best particle in the neighborhood; (a) Select the particle with the best fitness from the neighborhood of each particle as the best particle in the neighborhood; (b) When updating the particle position, a new variable G is introduced, which represents the influence of another neighboring particle that the current particle may be affected by. The calculation method of G is as follows: randomly select a different particle as a candidate solution of G. If the individual optimal fitness value of the particle is better than the individual optimal fitness value of the current particle, then the individual optimal position of the particle is assigned to G; otherwise, the global optimal position is assigned to G. (c) Introducing the best particle in the neighborhood by multi-point crossover to obtain the local attractor; specifically: calculating the value of each dimension of the local attractor through a loop; in each loop, first generating a random number; according to the size of the random number, randomly selecting the individual best position, the new variable G or the corresponding dimension value of the best particle in the neighborhood as the value of the corresponding dimension of the local attractor with equal probability; if the random number is less than 1 / 3, selecting the best position of the individual; if the random number is between 1 / 3 and 2 / 3, selecting the new variable G; if the random number is greater than or equal to 2 / 3, selecting the best particle in the neighborhood; after the loop is completed, returning a local attractor composed of all the local attractor dimension values; (d) Using the BQPSO algorithm, the probability value is calculated for each dimension; the following operations are performed on each bit of the substring corresponding to each decision variable of the local attractor: if the calculated probability value is greater than the random number, the bit is inverted; that is, if the bit was originally 1, it is set to 0; vice versa; and the substring is assigned to the corresponding substring of the updated particle; (e) After completion, all updated particle substrings are combined to obtain the new position; (f) Ensure that the updated particle positions satisfy specific constraints to maintain the validity of the problem solution space; (g) Use the FBA method to update the fitness value of the particle, and update the individual optimal solution and the global optimal solution; Step 6, iterative optimization: repeat steps 3 to 5 until the preset maximum number of iterations is reached to find the best strain transformation strategy; Step 7: Output optimization results: Output the global best particles as the optimal strain transformation strategy to guide strain transformation in the actual production process.

2. The metabolite optimization method based on the binary quantum particle swarm optimization algorithm and the flux balance analysis hybrid algorithm according to claim 1, characterized in that: Preprocessing is required before initializing parameters. The preprocessing steps are as follows: first, remove basic reactions required for growth from the list of potential knockout reactions, then exclude reactions that are not related to genes, and finally, do not consider biomass, ATP maintenance, and target reactions in the search space.

Images