Ind. Eng. Chem. Res. 1998, 37, 4709-4714
4709
Analysis and Optimization of Biochemical Process Reaction Pathways. 2. Optimal Selection of Reaction Steps for Modification Rau ´ l Conejeros and Vassilios S. Vassiliadis* Department of Chemical Engineering, Cambridge University, Pembroke Street, Cambridge CB2 3RA, U.K.
This is the second part of the study in the identification and manipulation of biochemical reaction limiting steps. The first part (Analysis and Optimization of Biochemical Process Reaction Pathways. 1. Pathway Sensitivities and Identification of Limiting Steps) establishes a link between biochemical reactions and the macroscopic fermentation process by using a set of amplification/attenuation factors which influence individually each reaction step. A sensitivity analysis uses these factors to establish the influence of each step on the process objective function for any processing mode. In this paper the optimal selection of the amplification/attenuation factors is examined, in formulations resulting in nonlinear programming. 1. Introduction This paper is the second part of the study in the identification and manipulation of biochemical reaction limiting steps. The first part established a systematic framework for identifying the limiting steps around the solution of an optimization problem, involving the performance of any fermentation process, focusing primarily in the dynamic processing mode. The framework is based on the evaluation of sensitivity measures of the overall kinetic rates over the objective function at the final process time. The current part of this work deals with the optimal selection of the amplification/attenuation factors of a subset of the enzymatic steps. This set is comprised of the steps exhibiting the largest sensitivity as identified with the methodology of the first paper. Of course, as these sensitivities are derived locally about an optimal solution, steps rejected from consideration may become influential at some distant operating point. Nonetheless, we may expect that it is desired to base information around a local point as large excursions may lead to totally unpredictable behavior of the kinetics in real situations. In the literature a number of approaches have been proposed on the identification of particular enzymatic steps and how to influence them in order to effect an overall improvement in the performance of a microorganism. The flux analysis approach has been used to estimate the magnitude of the reaction rates and variation of the concentration of intracellular and extracellular metabolites. Measurements are taken at some stages of the fermentation process so that a set of controlling rates can be determined.1,2 Optimization of the productivity of the reaction network has been done3-5 so that the maximum capacities of production of a pathway can be determined. These approaches require kinetics of the microbial system. Optimal modification of pathways6,7 was introduced by considering the influence of metabolites and enzymes over the reaction rates. This constitutes a new methodology in the improvements of the global productivity and modification of a biochemical reaction net* Author to whom correspondence should be addressed. E-mail:
[email protected].
work; the most significant consequence of their work is that a relaxation of the network structure is introduced in terms of the modification of its regulations. All of these approaches are set on a pathway modification in a pseudo-steady-state situation,8-12 assuming a continuous flow of material from substrate to product. None of these optimization approaches has considered dynamic aspects in the modification of pathways. Linearization of the system, according to the loglinear approach presented by Savageau,13-15 has been extensively used because of the high nonlinearity of the original kinetic expressions. Hatzimanikatis et al.16 have shown that there are only small deviations using this approximation. Therefore, the mixed-integer linear programming (MILP) problems solved for network modifications at steady state by Hatzimanikatis et al.,6,7 using the log-linear version of the kinetics, are likely not to be severely influenced as well by this approximation. The uncertainty of information within pathway models is significant, and that is why simplified models are usually adopted. Considering modification of a single step of the pathway is unlikely to have a big impact on the improvement of productivity. The uncertainty of information, or indeed its lack, prohibits the proposal for many modifications of steps simultaneously. Therefore, there is a tradeoff in optimization studies, where one or two steps may be unlikely to have an effect while too many may be either unreliable or impossible to implement. The current work is centered in the identification and modification of the limiting reaction steps within dynamic fermentation processes. The key differences proposed are as follows: 1. Consideration of amplification/attenuation factors for the key steps using full kinetics. 2. Delegation of the details of how to modify individual activation/inhibition of enzymes to a lower level analysis. 3. Modification of a small subset of reactions within an overall pathway. 4. Consideration of the process and the reaction network together. 5. Estimation of the minimal change in the pathway for a given productivity target.
10.1021/ie980411c CCC: $15.00 © 1998 American Chemical Society Published on Web 11/12/1998
4710 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
2. General Optimization of Amplification/ Attenuation Factors In the optimization problems of the first paper (OCP2), we can now consider releasing some of the parameters θi and treating them as optimization variables. The most influential steps can be selected according to
θi: i ∈ Sθ
(1)
where Sθ is the set of variables θi which induce the largest absolute sensitivity to the objective function; i.e.
| | ∂C* ∂θi
>
i∈Sθ
| | ∂C* ∂θl
(2) l∉Sθ
The dimension of the set Sθ is a user specification. This choice is then
dim[Sθ] ) Nθ e NR
(3)
The optimization of the amplification/attenuation factors involves the optimal selection of all θi belonging to set Sθ. Accordingly, the following bound constraints are appended to the original formulation of OCP-2, making all θi, i ∈ Sθ, additional design variables of the problem:
θL e θi e θU, i ∈ Sθ θi ) 1, i ∉ Sθ
(4)
case, where θi ) 1, the following objective can be minimized:
The bounds on θi are such that
0 < θL < 1 U
θ >1
Figure 1. Metabolic pathway optimization scheme.
(5a)
NR
min p
(5b)
Figure 1 shows the schematics of the metabolic pathway step optimization, which is set over the optimal control envelope proposed in the first part of this work. The restriction of the set of enzymes to modify Sθ may still contain many enzymes so that a restriction to a number of them for simultaneous modification should be imposed, e.g., modify only up to three enzymes, in which case the problem is combinatorial and a mixedinteger nonlinear programming (MINLP) problem should be solved. In the example case studies, it is demonstrated that using the ranking of enzymatic steps, provided by the sensitivity preprocessing stage, it is possible to generate very good optimal solutions approaching the fully relaxed case (all of the θi are free for optimization), through combination of a few steps at a given time, using the sensitivity ranked list. 3. Optimization-Based Ranking of Enzymatic Steps An alternative to ranking the steps, using their sensitivities, is by considering an optimization problem minimizing the modification of step amplification/ attenuation factors, in order to achieve a given productivity target. Thus, to increase the production of a product PROD of an unmodified fermentation process
(θi - 1)2 ∑ i)1
(6)
subject to the constraint
MPROD g R(MPROD)θ/i)1
(7)
(and subject to the dynamic process model); NR is the number of reaction steps of the pathway, and R > 1 is the desired increase factor in the the optimal amount / . Paof PROD produced in the base case (M)PRODθ i)1 rameter p is a sufficiently large penalty factor to ensure that the target constraint (eq 7) is satisfied with minimum modification of the amplification factors from the value of 1. A small value of R, i.e., slightly greater but close to unity, will highlight the direction in which all of the θi should move to achieve an increase in productivity. Thus, an immediate ranking of how much influence each step has is provided. For cases where the number of reactions is indeed very large, then optimization of many θi simultaneously might be too difficult computationally with existing optimization solvers. 4. Example Case Studies The model for the case studies in this paper is the same as that presented in Appendix B of the first part of this paper (Analysis and Optimization of Biochemical Process Reaction Pathways. 1. Pathway Sensitivities
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4711
and Identification of Limiting Steps). The model considers a part of the biochemical reaction pathway of Penicillium chrysogenum17 and matches the experimental data presented by Jørgensen et al.2,18,19 The setup in part 1 of this series is such that the process is considered as a fed-batch dynamic fermentation. The task is to maximize the production of penicillin V. In the first part of this work the sensitivity analysis was carried out, showing that the relaxation of the bottlenecks of the pathway can improve the performance of the fermentation. Because of the limitation in the improvement of the productivity with the amplification of a single step, a set of reaction steps must be considered for simultaneous modifications. The probability of a negative effect in the cellular viability can increase with the modification of a large number of steps in the pathway. To this end, the table of steps ranked according to their sensitivity magnitude can be used to select only a limited number of them for modification as indicated in section 2 of this paper. A number of case studies are constructed in the following subsections. All computational results were obtained using the equation oriented simulator gPROMS and the related general optimization package gOPT20 on a Sun SPARCstation 10 computer. 4.1. Optimization of θi in Order To Maximize Productivity. The objective function is set as the maximization of the amount of penicillin V produced by choosing an optimal control feed rate and modifying all, or a suitably chosen subset, of the enzymatic step amplification/attenuation factors θi. An optimization base case was run for two different final fermentation times, 110 and 220 h, where all of the values of θi were set to 1. The corresponding optimal penicillin V production levels obtained were 1142.16 and 2402.86 mmol for 110 and 220 h, respectively. The system was then optimized for a final time of 220 h. In this run the amplification/attenuation factors θi were all allowed to vary within lower and upper bounds set to 0.4 and 2.5, respectively. The result was that all θi for steps 1-5 were set to the upper bound of 2.5, while all steps 6-10 were taken to their lower bound of 0.4. The optimal productivity found for this case is 6151.03 mmol of penicillin V. Following this case study, we turn our attention to identifying a small subset of reactions for modification, instead of modifying every step in the pathway as in the previous case. If we choose to modify just three reaction steps for optimization, for nine reaction steps we can have 84 possible combinations to optimize. At this point we can predict that the best results of the optimization will be obtained for the combination of reaction steps with the higher sensitivities. This is a much simpler alternative to the use of a MINLP approach, which in our formulation with dynamics would be a mixed-integer optimal control problem (MIOCP), a far more complex problem to tackle with existing optimization techniques. Table 1 shows the result of the optimization for some of the triplets for a final fermentation time of 220 h. The results can be compared to a base case where all θi were allowed to vary within a range from 0.4 to 2.5. The results of the groups of reactions agree with the ranking of the reaction steps according to their sensitivities (Table 3 of part 1 of this work). In fact, the triplets were generated by including in the subset of three reactions
Table 1. Maximal Production of Penicillin V for Some Modified Reaction Step Triplets of the Pathway for a Final Time of 220 h triplet
optimal θi
MPenV (mmol)
base production (%)
all θ ) 1 objective (%)
base (2, 1, 5) (2, 1, 4) (2, 1, 8) (1, 5, 4) (1, 5, 8) (1, 5, 9)
(2.5, 2.5, 2.5) (2.5, 2.5, 2.5) (2.5, 2.5, 0.4) (2.5, 2.5, 2.5) (2.5, 2.5, 0.4) (2.5, 2.5, 0.4)
6151.03 5983.91 3073.53 2647.96 2541.35 2588.58 2516.57
100.00 97.28 45.97 43.05 41.32 42.08 40.91
255.99 249.03 127.91 110.20 105.76 107.73 104.73
the next best step to the worst one, going down on the sensitivity rank table. As can be seen in Table 1, through this procedure the most important triplet was chosen and its subsequent optimization produces 97.28% of the optimal amount of penicillin V of the case with all θi free. For the next triplets the productivity drops to ∼45% of the case with all θi free. In the case of triplet (1, 5, 4), the productivity obtained is below the value of the next triplets in the list (Table 1) produced by the sensitivity ranking procedure. This does not agree with the expected optimal production based on the ranking. The procedure does not assure global optimal points, and the results are local minima of the corresponding MIOCP allowing only three modifications of the nine reaction steps. Nonetheless, the discovery of the first triplet in our case study indicates that the procedure is a very useful tool to focus the pathway modification work to a few important reactions. 4.2. Minimization of Modifications To Achieve a Given Productivity Target. The sensitivity ranking of reactions was shown in the first part of this paper as a way to rank reactions according to their importance on the performance index, and they were used in the previous subsection to identify subsets of reactions for modification. In this section we use the idea of minimizing the deviation of the θi factors of each reaction in order to achieve a fixed productivity lower bound. An optimization problem can be set up to rank the greatest deviation from θi ) 1, with an objective of minimizing the overall modifications of the amplification/attenuation factors subject to a fixed lower bound on the yield. The new optimization problem considered is as follows: NR
(θi - 1)2 ∑ i)1
min 100
(8)
subject to the constraint
MPenV g R(MPenV)θ* i)1
(9)
(and subject to the dynamic process model); NR is the number of reaction steps of the pathway, MPenV is the product amount, and R is the desired increase factor in the optimal amount of penicillin V produced over the base case value M/PenV. Two runs were set up: (i) Minimal changes in θi for tf ) 110 hachieving at least the same production level as the base case run for tf ) 220 h, i.e., 2402.86 mmol of penicillin V. (ii) Minimal changes in θi for tf ) 220 h achieving twice the production level of the base case run for tf ) 220 h, i.e., 2 × 2402.86 ) 4805.72 mmol of penicillin V. As can be seen in Table 2, the main changes are in the forward reaction steps 1, 2, and 5. Reaction step 4
4712 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 Table 2. Minimal Deviation of θi for All of Them Free for the Two Cases at tf ) 110 and 220 h: List of Optimal θi and Their Deviations from 1
Table 5. Sensitivities When Just Reaction Steps 1, 2, and 5 Are Modified for a Final Time of 220 h: Maximal Productivity and Minimal Deviation Cases
final time (h)
sensitivities (mmol)
110
220
i
θi
θi - 1
θi
θi - 1
1 2 3 4 5 6 8 9 10
2.017 28 2.119 51 0.985 292 1.179 10 2.173 09 0.998 221 0.926 671 0.991 385 0.993 460
1.017 28 1.119 51 -0.014 71 0.179 10 1.173 09 -0.001 78 -0.073 33 -0.008 61 -0.006 54
1.768 76 2.073 20 0.984 219 1.154 13 1.958 08 0.981 267 0.939 044 0.993 258 0.996 561
0.768 76 1.073 20 -0.015 78 0.154 13 0.958 08 -0.018 73 -0.060 96 -0.006 74 -0.003 44
Table 3. Minimal Deviation of θi for Reaction Steps 1, 2, and 5 for the Two Cases at tf ) 110 and 220 h: List of Optimal θi and Their Deviations from 1
θi bounds (0.4-2.5)
minimized θi change
1 2 3 4 5 6 8 9 10
80 608 -50 334 1614 0 -132 -2 -1
416 611 -45 118 703 -35 -95 -27 -23
Table 6. Ranking of Significance of the Pathway Reaction Steps According to the Absolute Normalized Sensitivity of Penicillin V Production in Optimized Fed-Batch Culture at a Final Time of 220 h reaction step
final time (h) 110
i
220
i
θi
θi - 1
θi
θi - 1
1 2 5
2.013 76 2.119 15 2.212 41
1.013 76 1.119 15 1.212 41
1.789 41 2.098 27 2.018 45
0.789 41 1.098 27 1.018 45
Table 4. Sensitivities for All θi Free at a Final Time of 220 h: Maximal Productivity and Minimal Deviation Cases sensitivities (mmol) i
θi bounds (0.4-2.5)
minimized θi change
1 2 3 4 5 6 8 9 10
510 1868.8 1.6 -0.2 20.4 -0.8 -0.6 -11.4 -5
617 882 -36 109 754 -42 -75 -28 -25
has been amplified, directing the system to the production of penicillin V, but the change is not as significant as that in steps 1, 2, and 5. The remaining reaction steps have changes that are not as significant. As expected, the changes suggested by the optimization reduce the rates of the lateral reactions. In most of the steps the changes are minimum so that the efforts to modify the pathway should be directed just to a few reactions. The next optimization problem solved involved the same problem as that described in this section, i.e., minimal modification for the same case studies described in items i and ii, restricting changes only to the most important steps identified were 1, 2, and 5. The results of these two cases are summarized in Table 3. 4.3. Sensitivity Analysis of Various Case Studies and Other Related Results. The last point to check for the optimization runs is the change in the sensitivities due to the modification of the reaction steps. This is an important evaluation as it will indicate how the order of significance of the steps has shifted following the proposed optimal modification. The sensitivities are summarized in Table 4 for the cases where all θi are left as optimization variables free to vary within bounds 0.4-2.5 and for the fixed final time tf ) 220 h. The first column shows the sensitivities for the maximization of penicillin yield. The second column shows the sensitivities of the optimization run
θ (1, 2, 5) base all θ free all θ free θ (1, 2, 5) free case bounds minimized free bounds minimized rank all θ ) 1 (0.4-2.5) change (0.4-2.5) change 1 2 3 4 5 6 7 8 9
2 1 5 4 8 9 3, 10 6
2 1 5 9 10 3 6 8 4
2 5 1 4 8 6 3 9 10
5 2 4 8 1 3 9 10 6
5 2 1 4 8 3 6 9 10
minimizing the change of the θi to double the productivity of the base case with all θi fixed to 1. The same procedure, for the same two objective functions, but for modifications allowed only for reaction steps 1, 2, and 5 is shown in Table 5. The first column shows again the sensitivities for the reaction steps, where 1, 2, and 5 were allowed to vary between 0.4 to 2.5. The second column shows the deviations of the same reactions steps for the restricted modification of θi situation. All of these results show clearly a change in the value of the sensitivities, as expected, and the order of importance of the reaction steps. The change in the importance is due to the saturation of the sensitivity of some of the steps so that they are not as sensitive as those in the base case, where all θ ) 1. A rank of the reaction steps can be seen in Table 6 for all of these runs. Finally, it is interesting to provide an indication of the optimal control profiles for the feed rate of the optimization cases. Some of these profiles are shown in Figure 2, where it can be seen that they bare great resemblance with each other. 5. Conclusions and Future Work The first part of this work focuses on the derivation of biochemical pathway sensitivities and how they can be used to highlight the most important reaction steps for a particular production task. This part complements the sensitivity framework by using a smaller set of steps to consider for manipulation. The sensitivities, although local around an optimal operating point for the process, are used to identify reaction step bottlenecks and restrict proposed modifications to only the most important steps.
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4713
Figure 2. Feed-rate profiles for optimized reaction steps: (a) All θi allowed to vary between 0.4 and 2.5. (b) All θi deviations minimized to double productivity of the base case. (c) θi for steps 1, 2, and 5 allowed to vary between 0.4 and 2.5. (d) Minimized change of θi for steps 1, 2, and 5 to double productivity of the base case.
Obtaining sensitivities for several final processing times enables consideration of the rate of production as well as the final value of the objective function. Following this procedure, the most important reaction steps can be identified as those inducing the largest sensitivity to the objective in absolute value, and a small set of these can be selected for genetic engineering manipulation. The framework proposed here is a tool that maintains an overview of the pathway within any type of process and kinetics. By going through the levels of the hierarchy (sensitivities, step ranking, and modification of steps), one can extract valuable information for a pathway at the best possible set of operating conditions. In this respect, this work is original as it proves the impact of each reaction step on the optimal objective function value, including process dynamics and nonidealities. The hierarchical approach effectively reduces the complexity of the problem by increasing the resolution in the nested subtasks. Without losing generality, it can point to a subset of reactions that may be manipulated simultaneously for the best improvement in process performance. The issue of enzymatic activation/ inhibition by metabolites could be viewed as the final level in this hierarchy. This was not considered in this work. Using information based on the process and as much detailed kinetics as possible is a far better qualified
guess than arbitrary manipulations. This approach, coupled with expertise in biochemical reactions, can aid the genetic engineering task by reducing the number of trial manipulations required to improve the performance of a microorganism within general fermentation processes. Nomenclature p ) penalty factor MPenV ) cumulative mass of penicillin V (mmol) MPROD ) cumulative mass of product (mmol) NR ) no. of reaction steps in the pathway Nθ ) no. of amplification/attenuation factors PenV ) penicillin V (mM) Sθ ) set of amplification/attenuation factors θi R ) desired increase factor in the optimal amount of production θi ) attenuation/amplification factor for each reaction step i tf ) total fermentation time (h) Superscripts L ) lower bound U ) upper bound
Literature Cited (1) Delgado, J.; Liao, J. C. Inverse flux analysis for reduction of acetate excretion in Escherichia coli. Biotechnol. Prog. 1997, 13, 361.
4714 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 (2) Jørgensen, H.; Nielsen, J.; Villadsen, J. Metabolic flux distributions in Penicilliun chrysogenum during fed-batch cultivations. Biotechnol. Bioeng. 1995a, 46, 117. (3) Reagan, L.; Bogle, I. D. L.; Dunnill, P. Simulation and optimization of metabolic pathways. Comput. Chem. Eng. 1992, 16S, S237. (4) Reagan, L.; Bogle, I. D. L.; Dunnill, P. Simulation and optimization of metabolic pathways. Comput. Chem. Eng. 1993, 17, 627. (5) Voit, E. O. Optimization in integrated biochemical systems. Biotechnol. Bioeng. 1992, 40, 572. (6) Hatzimanikatis, V.; Floudas, C. A.; Bailey, J. E. Optimization of regulatory architectures in metabolic reaction networks. Biotechnol. Bioeng. 1996a, 52, 485. (7) Hatzimanikatis, V.; Floudas, C. A.; Bailey, J. E. Analysis and design of metabolic reaction networks via mixed-integer linear optimization. AIChE J. 1996b, 42, 1277. (8) Galazzo, J. L.; Bailey, J. E. Fermentation pathway kinetics and metabolic flux control in suspended and immobilized Saccharomyces cerevisiae. Enzyme Microb. Technol. 1990, 12, 162. (9) Galazzo, J. L.; Bailey, J. E. Errata. Enzyme Microb. Technol. 1991, 13, 363. (10) Schlosser, P. M.; Holcomb, T.; Bailey, J. E. Determining metabolic sensitivity coefficients directly from experimental data. Biotechnol. Bioeng. 1993, 41, 1027. (11) Schlosser, P. M.; Riedy, G.; Bailey, J. E. Ethanol production in baker’s yeast: Application of experimental perturbation techniques for model developement and resultant changes in flux control analysis. Biotechnol. Prog. 1994, 10, 141. (12) Torres, N. V.; Voit, E. O.; Glez-Alco´n, C.; Rodriguez, F. An indirect optimization method for biochemical systems: Description of method and application to the maximization of the rate of ethanol, glycerol, and carbohydrate production in Saccharomyces cerevisiae. Biotechnol. Bioeng. 1997, 55, 759.
(13) Savageau, M. Biochemical system analysis. i. some mathematical properties of the rate law for the component enzymatic reactions. J. Theor. Biol. 1969a, 25, 365. (14) Savageau, M. Biochemical system analysis. iii. dynamic solutions using a power-law approximation. J. Theor. Biol. 1970, 26, 215. (15) Savageau, M. Biochemical system analysis. ii. the steadystate solutions for an n-pool system using a power-law approximation. J. Theor. Biol. 1969b, 25, 370. (16) Hatzimanikatis, V.; Bailey, J. E. Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models. Biotechnol. Bioeng. 1997, 54, 91. (17) Pissara, P. N.; Nielsen, J.; Bazin, M. J. Pathway kinetics and metabolic control analysis of a high-yielding strain of Penicillium chrysogenum during fed-batch cultivations. Biotechnol. Bioeng. 1996, 51, 168. (18) Jørgensen, H.; Nielsen, J.; Villadsen, J.; Møllgaard, H. Analysis of penicillin v biosynthesis during fed-batch cultivations with a high-yielding strain of Penicillium chrysogenum. Appl. Microbiol. Biotechnol. 1995b, 43, 123. (19) Nielsen, J.; Jørgensen, H. Metabolic control analysis of the penicillin biosynthetic pathway in a high-yielding strain of Penicillium chrysogenum. Biotechnol. Prog. 1995, 11, 299. (20) gPROMS User’s Guide, Version 0.1; Process Systems Enterprise Ltd.: London, 1996.
Received for review June 24, 1998 Revised manuscript received September 9, 1998 Accepted September 10, 1998 IE980411C
