Ind. Eng. Chem. Res. 1998, 37, 4699-4708
4699
Analysis and Optimization of Biochemical Process Reaction Pathways. 1. Pathway Sensitivities and Identification of Limiting Steps Rau ´ l Conejeros and Vassilios S. Vassiliadis* Department of Chemical Engineering, Cambridge University, Pembroke Street, Cambridge CB2 3RA, U.K.
This work focuses on establishing 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 based on these factors yields a novel measure of reaction step influence on the overall process objective. These measures are generalized to any type of processing mode, primarily focusing on dynamic processes, for which it is shown how the sensitivities can be used to identify reaction network bottlenecks. 1. Introduction
r)
The knowledge of the metabolic routes can help in the determination of the cellular productive capabilities and ways of improving them. Thus, the maximum theoretical yield can be estimated from the forward chemical reactions from substrate to product. In a cellular system this estimated value is reduced due to the multitude of intermediate steps and lateral reactions in the pathway.1,2 If some of these reactions could be avoided or slowed down, the product yield could be improved. Because of the expansion of knowledge of genetic sequences and transcription controls, it is possible to modify the production of metabolites and regulation of the control system of a cell. To this end, mathematical tools have been developed to estimate the consequences of the manipulation of the pathway and in this way to analyze at the same time the viability and stability of the expected result.3-7 In the current literature a number of approaches have been proposed for the analysis of biochemical reaction pathway networks. These can be categorized under rate flux analysis studies, a combination of flux analysis and optimization, and recently the use of optimization methods exclusively, such as linear programming (LP), nonlinear programming (NLP), and mixed-integer linear programming (MILP). Most of the approaches have considered the reaction kinetics for a partial set of the complete pool of cellular reactions, trying to optimize the amount of metabolites inside the cell around the steady state of the system. The metabolic flux analysis is based on a stoichiometric model that can be applied to cell reactions. The metabolic flux can be defined as the molar rate of production of a metabolite. The stoichiometry for all of the cellular reactions can be specified as in8
A‚X ) O
(1)
where the matrix A contains all of the stoichiometric coefficients and the vector X contains all of the metabolites involved in the metabolic pathways. Considering a forward rate vector v, the net production rates for the intracellular metabolites can be expressed as * Author to whom correspondence should be addressed. E-mail:
[email protected].
dX ) ATv - µX dt
(2)
The first term of this equation gives the consumption and production of metabolite i, and the second is a dilution term due to an increase of the total amount of biomass during the fermentation process. Thus, the forward rate can be calculated as
v ) (AT)-1(r + µX)
(3)
Usual steady-state analysis9 considers a continuous flow of material through the pathway, and the center of the analysis is the set of reactions. The dilution term and consumption of species in growth are not usually considered. In the conventional flux analysis, the variation of extracellular species such as substrates, products, and gases is measured during the fermentation process. Considering the biochemical pathway and stoichiometry, the internal fluxes can be calculated using the steady-state approach. Depending on the regulatory controls inside the cell, the flux values are not constant during the fermentation time10 so that the steady-state assumption is not always valid in all the stages of fermentation. The available models for steady-state situations11-13 do not describe all of the extent of a dynamic fermentation process.14 Recent publications have considered the optimization of metabolic pathways3-7,9,15 in order to obtain an optimal network modification for the overproduction of a specific metabolite. These papers are based on the linearization proposed by Savageau,16-19 starting from a power law description of the complex nonlinear reaction kinetics. The main consideration in the optimization is the modification of the regulations over some of the steps of the pathway by the substrates, products, or any other metabolites involved. The results of this approach provide an indication of the necessity for overproduction or underproduction of an enzyme, or the modification of the structure of the enzyme involved, in order to change the influence of metabolite concentrations. A recent approach9 has considered the effect of the amplification or attenuation of some of the steps in the reaction network, in order to obtain the optimal over-
10.1021/ie980410k CCC: $15.00 © 1998 American Chemical Society Published on Web 11/12/1998
4700 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
expression or underexpression of an enzyme, giving a target for a possible genetic engineering manipulation. The stability of the optimized reaction network is studied by Voit3 in order to assure the consistency of the proposed change through a transient response of the system to a step perturbation. Hatzimanikatis et al.20 made a comparison between the linearized power law model and the nonlinear Monod kinetic-based model, showing that there is no significant deviation in the linearized system for perturbations in the network of reactions. The work of Hatzimanikatis et al.6 is a new approach in the area of model-guided genetic engineering predictions. Using mixed-integer linear programming, it identifies points where the network could benefit by genetic manipulations so given enzymes may become more active by overexpression or inhibitory/activation loops are modified (or introduced) artificially so that the given objective is maximized. Mauch et al.21 have recently proposed the study of the reaction network as a part of the metabolic system, considering the whole fermentation process in the analysis of a general framework. Their focus was steady-state operation. The issue of dynamics is still open, and no work has actually considered their implications on the proposed changes in the underlying biochemical reaction network. The classical metabolic control theory is performed over a particular point of the fermentation, assuming that the production and consumption of metabolites involved are balanced.7,22,23 The perturbations of any state property or parameter give an idea of the effect of its variation over the system behavior. This effect on the material flux J by the enzyme variation E is known as control coefficient analysis.
C)
∂J E ∂E J
(4)
In an analog way the influence of the metabolites over the reaction rate is known as elasticities analysis. Control coefficients and elasticities can be estimated from experimental data or calculated directly from kinetic expressions, when kinetic information is available. Once obtained, the analysis of the system turns in a pure linear structural analysis around the local region of the pseudo steady state.22,23 Improvements suggested by this method are merely local estimations, and the control coefficients can change in different stages of the fermentation, so depending on the information with which the coefficients were calculated, the result may not be valid for the whole process. We propose a more extensive approach considering the whole fermentation process and the effect of the changes on the metabolic pathway over the overall productivity of the system. Its result should be global for the fermentation process, considering all of the possible changes of the control coefficients. The need for kinetics might be viewed as a restriction of the applicability of the proposed methodology. Nonetheless, the question of whether kinetics should be available does not weaken the theoretical contribution and potential of the technique. In this work we present an automated procedure to estimate the productivity expectations in a fermentation system through an analysis of the metabolic reaction network. The sensitivity of the system to a perturbation
of every reaction step will provide the location of the bottlenecks of the system. This may serve as a guide to a possible genetic modification of a biochemical pathway, in order to improve its productivity characteristics. In addition, the procedure gives an estimation of the expected improvement, as is demonstrated in the companion paper (Analysis and Optimization of Biochemical Process reaction Pathways. 2. Optimal Selection of Reaction Steps for Modification). 2. Dynamic Model and Optimal Control Problem The general description of the fermentation process will give us the framework to determine the capabilities of a biochemical pathway for maximum productivity for a fixed processing time. The achievement of this objective is limited by the restrictions imposed by the structure and regulations of the specific metabolic pathway, which generate the unique dynamics and behavior for the microorganism. The control variables provide modifications of some of the external conditions and can be used to drive the fermentation process to maximize the desired objective. The optimization task results in an optimal control problem (OCP), which can be described according to the following expression:
OCP-1 max x(‚),r(‚),u(‚),v
(5a)
ri(t) ) fi(x(t),u(t),v), i ) 1, 2, ..., NR
(5b)
x˘ j(t) ) gj(x(t),r(t),u(t),v), j ) 1, 2, ..., NX
(5c)
subject to
x(t0) ) x0 t0 e t e tf
fixed initial conditions fixed final time
3. Sensitivity Framework In designing this framework, we set the following requirements: 1. Automatically identify limiting reaction steps by introducing a measure which can be evaluated for all of the steps in any reaction pathway. 2. The sensitivity measures introduced should be applicable for entire processes, treating the biochemical pathway as any other reaction engineering system, and maintaining an overview of the outer process control and design parameters. 3. Be able to look at any type of processing mode, i.e., batch and fed-batch operation, which are dynamic in nature and involve time-varying feed rates; the latter calls for an optimal control formulation. 4. Use kinetics of arbitrary complexity and nonlinearity; in principle, we seek a methodology that is able to look at a very large network in its original kinetic form and identify limiting steps. 5. The sensitivity formulation should not only be able to identify definitively the limiting steps but also be used consistently in a predictive mode to select a number of steps which should be considered for simultaneous genetic engineering manipulation. Mauch et al.21 have proposed a similar framework to analyze the metabolic system and to show in a general example the dependence of the metabolic output with
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4701
the sensitivity of the network in a steady-state condition. The difference in our work is the determination of the sensitivity of the optimal production of the system viewed as a dynamic process. The amplification or attenuation of steps, at any level of the production sequence, could have an effect on the substrate demand or on the internal regulation of the pathway, so that the optimal control profile will change completely. In order to determine the controlling bottlenecks in the optimized system, an artificial factor θi was applied to each step of the reaction pathway. The most important key features we require are the ability to look at reactions and processes in tandem, not in isolation, and be able to tackle very large-scale networks without resorting to kinetic expression simplifications. The sensitivity problem is formulated according to the next set of equations:
OCP-2 max C(x(tf)) x(‚),r(‚),u(‚),v
(6a)
subject to ri(t) ) wi(θi,x(t),u(t),v) ) θifi(x(t),u(t),v), i )1, 2, ..., NR (6b) x˘ j(t) ) gj(x(t),r(t),u(t),v), j ) 1, 2, ..., NX x(t0) ) x0
(6c)
fixed initial conditions
t0 e t e tf
fixed final time
where all variables appearing are as defined in eq 5a-c except for a new set of functions wi indicating the modified kinetics by the new parameters θi. The latter form an artificial parameter vector, of size equal to NR, and they are introduced as multiplication factors of the overall reaction rates for each individual step in the pathway. The original 7 pathway (and optimization problem) is obtained by setting θi ) 1, where i ) 1, 2, ..., NR. Similar to the analysis presented by Stephanopoulos et al.9 which is applicable to reaction pathway networks in the steady-state case, the sensitivity framework is built over the optimization envelope as shown in Figure 1. To study the effect of each individual rate expression on the objective function, we form the following augmented performance index (similar to the Lagrangian in steady-state cases) by adjoining all of the constraints to the objective and studying its perturbations at the optimal solution point:
J ) C*(tf) +
∫tt [λ*(t)]T(w*(t) - r*(t)) dt + ∫tt [µ/L(t)]T[g*(t) - x˘ (t)] dt f
0
f
0
(7)
where λ are Euler-Lagrange multipliers for the reaction rate equations (5b) and µL are Euler-Lagrange multipliers for the material balance equations (5c). The main result, proven in Appendix A, relates the sensitivity of the objective function to the individual parameters θi and is
|
∂C* (t ) ∂θi f
) at optimum
∫tt λ/i (t) fi(x*(t),u*(t),v*(t)) dt, f
0
i ) 1, 2, ..., NR (8)
Figure 1. Sensitivity analysis of the biochemical reaction pathway.
The importance of this result is that it extends the physical meaning of the Euler-Lagrange multiplier for the entire duration of the dynamic process and gives and an overall measure of the effect each individual reaction step has on the objective function through its rate at each time. It is evident that in the dynamic case the θ factors encapsulate very compactly the effect of each reaction step, while such a study would be extremely difficult by just looking at the profiles of the Euler-Lagrange multipliers. In fact, in many optimal control approaches,24-27 their profiles are not readily available and extra effort would be required to obtain them. It is finally noted that use of the sensitivity approach in chemical kinetics is a standard means of practice in order to classify the importance of various steps on the evolution of concentration profiles and yields as functions of time.28,29 The major difference in our approach is that (a) an artificial parameter is introduced which affects the reaction rates at all times, and later on is chosen optimally, and (b) the sensitivity is evaluated with respect to a process objective function at the final time. The parameters θi treated as amplification/attenuation factors can attain also a physical meaning by considering the derivation of enzymatic kinetics. These kinetics are assuming a first-order reaction of the enzyme and substrate and taking a fast equilibrium of the enzyme-substrate intermediate complex, according to the Michaelis-Menten hypothesis.30 The latter equilibrium reaction is proportional to the enzyme concentration; thus, the implication of the θi factors is a
4702 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 Table 1. Computational Statistics for Base Case Optimization Runs final time (h) total CPU time (s) optimizer statistics no. of optimization iterations no. of linesearch steps infeasible linesearch steps
Figure 2. Reaction pathway for penicillin V production in Penicillium chrysogenum.
requirement of scaling up or down the enzyme concentration by a fixed factor at all times. 4. Example Case Study This example is based on the model for Penicillium chrysogenum in fed-batch cultures.31 Figure 2 shows the schematic reaction pathway of this cell, which will be used in this work. The model equations for each reaction step of the pathway, as well as the dynamic material balance for each intracellular and extracellular metabolite, are included in Appendix B. The model describes the variation of the metabolite concentration involved in the penicillin-V production pathway along the fermentation time. It assumes that the concentration of the precursor amino acids is constant.31 This is a limitation of the model, because it does not connect completely the growth kinetics and overall fermentation material balance equations with the penicillin production. Still, our methodology is up to now one of the few approaches trying to describe the fermentation process and modeling the kinetic behavior of part of a metabolic pathway. Whenever a better kinetic description of the pathway is available, this approach is applicable. Following the methodology proposed to determine the sensitivities of the system, each reaction rate of the pathway was multiplied by a factor θi, where i ) 1, 2, ..., 10. The simulations involved were solved using the general process modeling system package gPROMS and the optimization runs were carried out using the general optimization package gOPT.32 The settings for the computational runs were according to the experimental conditions given by Jørgensen et al.10,33 in previous works, and results of the simulations were checked 9 with the results presented by Pissara et al.31 The conditions considered are initial volume V0 ) 4 L, final volume Vf ) 25 L, initial concentration of glucose in the fermenter so ) 3 g/L, initial biomass concentration in the fermenter x0 ) 4g/ L, feeding glucose concentration sF ) 450 g/L, and maximum fermentation time of tf ) 220 h, and the flow rate was allowed to vary within a range from 0 to 1 L/h. Finally, the feed flow rate was considered as the control variable of the system which to be optimized over the duration of the process. Optimization runs were made for four final times, 55, 110, 165, and 220 h. For each reaction step a variation
55
110
165
220
477.3
1325.0
1046.8
1138.5
39 67 0
47 80 0
31 59 0
36 61 1
for ∆θi of 1% was used to evaluate the sensitivities. A base case was run for each final time where all θi ) 1. The initial profile for the substrate flow for these cases was set at 0 for each one of the 20 piecewise constant control elements used in the optimization runs. The sensitivities δXi/δθi were thus calculated by a finite difference scheme, using the base case as the reference value. For each perturbed θi optimization run, the flow profile obtained in the base case was used as a starting profile, in this way minimizing the CPU time required to converge them. The optimization run statistics for the base cases are shown in Table 1, with CPU times reported on a Sun SPARCstation 10 computer. The control profiles obtained for each fermentation time in the base cases, where all θi ) 1, are shown in Figure 3. The summary of the sensitivity results is shown in Table 2. Clearly, it can be seen that reactions 1, 2, and 5 are the most important in the control of the production of penicillin. An important point is the change of the bottleneck after 110h from reaction 5 to reaction 2. This means that a decision to improve the productivity of the system by modifying the expression of enzyme 2 only because of results obtained in an analysis for short fermentation times could have no significance for longer times. This agrees with the idea that steady-state results could drive to wrong conclusions because of the dynamic nature of fermentation processes. Pissara et al.31 simulated the control coefficients for the fermentation process using this model, assuming a pseudo steady state for each point of the process. They consider this not to be a big error because of the fact that the characteristic time is greater than 100 h, more than 3 times the turnover of the LLD-ACV pool, which is considered the metabolite that varies the most. Their results show that the significant control coefficients are the ones for the enzymes IPNS (reaction step 2), ACVS (reaction step 1), and PA (reaction step 8) and a change in the rank of the magnitude of the control coefficient is observed. That means changes in the corresponding flux distribution, which agrees with the experimental results presented in literature,10,34 where the calculated control coefficients exhibit a similar profile. According to the “summation theorem” the sum of the control coefficients should be 1 for a pathway23
CACVS + CIPNS + CPA ) 1
(9)
In this case because of the fact that just three of the control coefficients are not zero or negative, the change in the magnitude of one of the control coefficients over the others means a shift in the limiting step in the control of the flux of the whole pathway. So, as expected, this is perfectly shown in the results obtained through the sensitivity analysis. The difference is that the reaction steps obtained as more important are the
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4703
Figure 3. Optimal substrate flow-rate profiles obtained for fed-batch cultures of P. chrysogenum at four different final times: (a) tf ) 55 h; (b) tf ) 110 h; (c) tf ) 165 h; (d) tf ) 220 h. Table 2. Sensitivity of Penicillin V Production (δMPenV/ δθi(g)) in Optimized Fed-Batch Culture at Different Final Times final time (h) i
55
110
165
220
1 2 3 4 5 6 8 9 10
27.5 28.1 -1.9 22.5 135.7 -0.1 -6.4 -0.3 -0.2
127 195 -10 111 584 -1 -43 -2 -1
373 1044 -1 49 277 -2 -19 -5 -2
411 1552 -2 47 265 -3 -20 -6 -2
ones affecting the overall productivity of the fermentation process and not just the ones affecting the pathway. Reaction steps 1 and 2 are shown as important in both analyses. For the four final times considered, 55, 110, 165, and 220 h, the first pathway steps are ranked using the absolute value of the normalized sensitivities shown in Figures 4-6, so that in a large reaction system, the most significant steps to modify can be identified easily. The normalization was done considering a value of zero for the base case and 1 for the largest positive sensitivity. In this case, we chose to consider as significant the steps with a relative sensitivity greater than 10% of the highest observed. The ranking of the reactions for the four final times is shown in Table 3.
Figure 4. Normalized sensitivity of reaction steps 1, 2, 4, and 5 for four different final fermentation times.
The possible improvements in penicillin production due to the modification of a single step of the pathway were studied considering the ranking of significance obtained for six levels of amplification, 1, 10, 25, 50, 100, and 200%, and for two final times, tf ) 110 h and tf ) 220 h. As can be seen in Figure 7 for the case of reaction 2, amplification of the expression of the enzyme from 25 to 200% yields no significant improvement in the amount of penicillin V produced for tf ) 220 h. The effect on the amplification of the step has a notable effect until the reaction step is not limiting the dynamic
4704 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
Figure 5. Normalized sensitivity of reaction steps 3 and 8 for four different final fermentation times.
Figure 6. Normalized sensitivity of reaction steps 6, 9, and 10 for four different final fermentation times. Table 3. Ranking of Significance of the Pathway Reaction Steps According to the Absolute Normalized Sensitivity of Penicillin V Production in Optimized Fed-Batch Culture at Different Final Times final time (h) rank
55
110
165
220
1 2 3 4 5 6 7 8 9
5 2 1 4 8 3 9 10 6
5 2 1 4 8 3 9 6, 10
2 1 5 4 8 9 6, 10 3
2 1 5 4 8 9 3, 10 6
response of the pathway, with the bottleneck moving to another reaction step. The same is observed for reaction steps 1 and 5 with a less marked effect over the productivity of penicillin V. A similar case is shown for reaction 5 in Figure 8 for tf ) 110 h. In this case reaction steps 1, 2, and 4 have less effect on the productivity of penicillin V, due to a shifting of the bottleneck. The logical next step, exploring the capabilities of the network, is the optimization of the system fixing a higher productivity target so that the required modifications in the pathway can be quantified. Modification of a large number of reaction steps must be avoided as it can have unknown consequences in the viability of the cellular system, making the high productive objective hard to reflect real situations. A proposal for the modification of groups of reactions in a metabolic system is given in the companion paper.
Figure 7. Penicillin V production under different amplifications of the pathway reactions 1, 2, and 5 for a final fermentation time of 220 h.
Figure 8. Penicillin V production under different amplifications of the pathway reactions 1, 2, 4, and 5 for a final fermentation time of 110 h.
4.1. Summary of Findings in the Case Study. 1. For the final times tested a change in the rank of importance of the reaction steps was obtained, so that the bottlenecks of the pathway move at different stages of the fermentation. The novelty of the approach is that complex dynamic systems can be studied. 2. The magnitude of the sensitivities rank each reaction step according to its importance in the pathway. The highest sensitivity values indicate bottlenecks of the pathway and can be used as indicators for very large pathways. 3. To improve the product output of the pathway, a single step can be amplified or attenuated. We chose to modify all the reaction steps with a sensitivity greater than 10% with respect to the largest one. It is shown that the productivity can be improved up to a limit in the magnitude of the amplification. Beyond that value, there is no improvement, at which point another reaction step restricts the pathway. 4. This approach can handle kinetics in their natural form, indicating immediately which single-step modification will give the best optimal result without linearization or any simplification of the kinetics and includes the overall process dynamics under consideration. 5. Conclusions and Future Work The formulation presented for dynamic operation bottlenecks in biochemical reaction pathways offers a unique analysis framework which is more powerful than simple simulations, and where kinetic expressions are
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4705
available, it is more general than flux analysis studies. Since the optimization of the operation of a given process is considered, the reaction pathway is treated within the overall process and for the optimal operating and design conditions identifies the bottlenecks attributable to reaction steps. The formulation considers the fact that the bottlenecks in a network not only may shift in value from one fixed final time to another, but also indicate that another step becomes dominant in the network at that production policy. Thus, the prediction of genetic changes can include results relating to the “processing speed” of a microorganism for a given type of processing mode. The derivation of the results assumes no particular restrictions on the kinetic expressions. It is also noted that the evaluation of the sensitivities does not load the optimization problem (OCP) with extra parameters. If the Euler-Lagrange multipliers are available, a sensitivity calculation with respect to the θ factors can be solved at the optimal set of control profiles u(t) and design parameter values v. Alternatively, as was done in this paper, the sensitivities may be calculated directly by re-optimizing the problem subject to finite difference perturbations of θj one at a time, using the base optimal solution as a starting point to speed up the new optimization problems. The values of the sensitivities exist for any reaction step. They can be used to answer questions as to where the bottlenecks are in large pathways and by userdriven ranking to identify the most important steps to be considered for genetic manipulation. Through the use of sensitivities, the focus of attention can be directed to a smaller part of a very large network of reactions, without loss of an overview of the overall network. Shifting from one objective function to a different one is also a very interesting prospective, as one optimization and sensitivity analysis will reveal the shift in relative importance of the different reaction steps, highlighted by their sensitivity measures’ magnitude and sign. Finally we do note that one possible limitation is the fact that the sensitivities are evaluated at one solution point. Thus, a decision to drop a step from consideration for modification is based on local information, and some solutions may be lost. Nonetheless, the preprocessing phase which is based on the sensitivities will be able to reduce the scope to a few reactions that in real practice are always inspected for correctness. These issues are examined further in the accompanying paper. Nomenclature f ) vector of forward reaction rates (mM h-1) g ) vector of material balance equations for each of the NX compounds (mM h-1) kACV ) secretion constant of intracellular ACV (h-1) ki ) first order reaction rate constant (h-1) kIPN ) secretion constant of intracellular IPN (h-1) m ) maintenance coefficient (g g-1 h-1) mM ) maximum maintenance coefficient (g g-1 h-1) r ) net production rates r ) vector of reaction rates (mM h-1) s ) substrate, extracellular glucose concentration (g L-1) s0 ) initial concentration of substrate in the fermenter (g L-1) sF ) substrate concentration in the feed stream (g L-1) t ) time (h)
t0 ) initial fermentation time (h) tf ) total fermentation time (h) u ) vector of control variables v ) vector of forward reaction rates (mM h-1); vector of design parameters w ) modified reaction rates by the amplification/attenuation factor θi x ) biomass concentration (g L-1) x0 ) concentration of biomass in the fermenter at t ) t0 (g L-1) A ) stoichiometric coefficients’ matrix AAA ) R-aminoadipic acid (mM) AAT ) acyl CoA to 6-APA transferase ACV ) L-R-aminoadipyl-L-cysteinyl-d-valine (mM) ACVOUT ) extracellular concentration of LLD-ACV (mM) ACVS ) L-R-aminoadipyl-L-cysteinyl-D-valine synthetase AT ) acyltransferase C ) cost function; flux control coefficient CoA ) coenzyme A CYS ) cysteine (mM) E ) enzyme concentration (mM) F ) feed flow rate (L h-1) GSH ) glutathione (mM) HPA ) 8-hydroxypenillic acid (mM) IAH ) isopenicillin N aminohydrolase IAT ) acetyl CoA to isopenicillin N acyl transferase IPN ) isopenicillin N (mM) IPNOUT ) extracellular concentration of IPN (mM) IPNS ) isopenicillin N synthetase J ) flux intracellular metabolites for a given reaction (mM h-1) J ) dual function of the optimal control problem KAAA ) saturation constant of ACVS for L-R-aminoadipate (mM) KACV ) inhibition constant of ACVS for ACV (mM) KCYS ) saturation constant of ACVS for L-cysteine (mM) Ki ) inhibition constant of IPNS for glutathione (mM) Km ) saturation constant of IPNS for ACV (mM) Km,6APA-Poa ) saturation constant of AAT for 6-APA with phenoxyacetyl CoA (mM) Km,IPN ) saturation constant of IAH for IPN (mM) Km,IPN-Poa ) saturation constant of IAT for IPN with phenoxyacetyl CoA (mM) Km,Poa ) saturation constant of AAT and IAT for phenoxyacetyl CoA (mM) KPenV ) saturation constant of PA for penicillin V (mM) Ks ) saturation constant of the specific growth rate with substrate (g L-1) Ksm ) saturation constant for the maintenance coefficient expression (g L-1) KVAL ) saturation constant of ACVS for L-valine (mM) V ) volume (L) V0 ) volume in the fermenter at t ) t0 (L) Vf ) volume in the fermenter at t ) tf (1) VAL ) valine (mM) IPN ) isopenicillin N (mM) 8HPA ) 8-hydroxypenillic acid (mM) 6APA ) 6-aminopenicillanic acid (mM) NR ) total number of reactions NX ) total number of metabolites O2 ) oxygen concentration (mM) OPC ) 6-oxopiperidine-2-carboxylic acid (mM) PA ) penicillin amidase PACVS ) linear decay constant for ACVS activity (mM-1 h-1) PAT ) linear decay constant for AT activity (mM-1 h-1) PINPS ) linear decay constant for IPNS activity (mM-1 h-1) PenV ) penicillin V (mM) PoaCoa ) phenoxyacetyl CoA (mM) X ) metabolites’ concentration vector (mM) XACVS ) ACVS activity (mM)
4706 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 XAT ) acyl transferase activity (mM) XIPNS ) IPNS activity (mM) Y ) substrate to biomass yield (g g-1) Yo ) maximum yield of substrate to biomass (g g-1)
Table 4. Model Parameters
Greek Letters λ ) Euler-Lagrange multipliers for the reaction rate equations µ ) specific growth rate (h-1) µL ) Euler-Lagrange multiplier for the concentration variation of metabolites µM ) maximum specific growth rate (h-1) Fc ) cell specific volume (mL g-1) θi ) attenuation/amplification factor for each reaction step i Superscripts / ) value of variable at the optimal point T ) transpose matrix
Appendix A We present here the sensitivity analysis and prove the final result for OCP-2 (eq 6a-c). The EulerLagrange multipliers λi, where i ) 1, 2, ..., NR, for each of the reaction rate equations, and µLj, where j ) 1, 2, ..., NX for each of the balance equations, are introduced to yield the following dual function (same as eq 7):
J ) C(tf) +
∫tt [λ(t)]T(w(t) - r(t)) dt + ∫tt [µL(t)]T(g(t) - x˘ (t)) dt
parameter
value
units
Km,IPN-Poa Km,Poa Km,IPN Km,6APA-Poa Ki Km KAAA KCYS KVAL KACV KPenV k1 k2 k3 k4 k5 k6 k8 kACV kIPN PACVS PIPNS PATPoaCoa GSH O2
0.023 0.006 4 0.0093 8.9 0.13 0.63 0.12 0.3 12.5 2 17.77 16.85 4.03 1.95 13.74 0.065 0.4 0.05 0.023 615 330 420 0.006 3 0.1157
mM mM mM mM mM mM mM mM mM mM mM h-1 h-1 h-1 h-1 h-1 h-1 h-1 h-1 h-1 mM-1 h-1 mM-1 h-1 mM-1 h-1 mM mM mM
condition for a stationary point. Since the initial condition x0 is fixed, we obtain
|
∂C - µTL(tf) ) 0 ∂x t)tf
f
0
f
0
(A.1) λT(t)
The property of the dual function above is that at the optimal solution it should equal the value of the objective function:
J * ) C*|at constrained optimum
-λT(t) + µTL(t)
(A.2)
The equation of the dynamic dual function is first modified by integrating by parts the last term to yield
J ) C(tf) +
∫t [λ(t)] (w(t) - r(t)) dt + tf
T
0
∫t [µL(t)] g(t) dt + [µL(t)] X(t)|t)t - [µL(t)] X(t)|t)t + ∫tt [µ˘ L(t)]TX(t) dt (A.3) tf
T
T
T
0
0
f
f
0
The necessary conditions of optimality for OCP-2 are obtained by considering feasible perturbations in δv and δu to obtain the corresponding perturbations in δx and δr (the reason for the integration by parts was to avoid having perturbations in δx˘ ):
δJ )
(∂C∂x δx))| - (µ δx)| + (µ δx)| + ∂w ∂w δx + δu + δv - δr] dt + ∫ λ [∂w ∂x ∂u ∂v ∂g ∂g ∂g δx + δu + + δr dt + ∫ µ [∂g ∂x ∂u ∂v ∂r ] t)tf
T L
t)tf
T L
t)t0
tf T
t0
tf T L
t0
∫tt µ˘ TL δx dt f
0
∂w ∂g + µTL(t) + µ˘ TL(t) ) 0 ∂x ∂x
(A.4)
To obtain the necessary conditions, the multiplier values of δx, δr, δu, δv are all set to zero. This will result in the condition δJ ) 0, which is the necessary
δJ * )
[({
∂g )0 ∂r
(A.5) (A.6)
(A.7)
λT(t)
∂w ∂g + µTL(t) )0 ∂u ∂u
(A.8)
λT(t)
∂w ∂g + µTL(t) )0 ∂v ∂v
(A.9)
} )
∫{
tf ∂C ∂x ∂w - µTL + δθ + t λT(t) 0 ∂x ∂θ t)tf ∂x ∂x ∂w µTL(t) + µ˘ TL(t) δθ dt + ∂x ∂θ tf ∂g ∂r δθ dt + -λT(t) +TµTL(t) t0 ∂r ∂θ tf T ∂w ∂g ∂u λ (t) + µTL(t) δθ dt + t0 ∂u ∂u ∂θ tf T ∂w λ (t) (A.10) + δθ dt t0 ∂θ at the optimal solution
∫{ ∫{ ∫
}
]|
} }
These, together with the original differential equations (6c) of OCP-2, form the set of necessary conditions of optimality. Conditions (A.6) and (A.7) are the adjoint DAE equations, while conditions (A.8) and (A.9) are equivalent to the Hamiltonian conditions of classical optimal control.35 By considering perturbations in δθ about an optimal solution for u/(t) and v/ and their effect on the dual function, we obtain By the necessary conditions in eqs A.5-A.9, all the terms in braces in eq A.10 become zero, and using the definition of functions w in eq 5b, and eq A.2
Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998 4707
for δJ *, leads to the result
|
∂C* ) ∂θ at constrained optimum
∫tt [λ*]T‚f(x*,u*,v*) dt f
0
(A.11)
which element wise is
|
∂C* ∂θi
)
at constrained optimum
∫
tf / λ ‚f (x*,u*,v*) t0 i i
dt,
i ) 1, 2, ..., NR (A.12)
Appendix B The description of the system used as a case study is described in Pissara et al.31 Equations B.1-B.10 give the reaction rates ri. The units are given in mmol L-1 Cell h-1. The specific volume is 2.4 (mL g-1 dw). It is noted that, in the reaction steps below, reactions 3 and 7 are taken to have the same rates. In sensitivity tables, only the modification of reaction rate 3 will be reported, as modification of reaction 7 will be forced to be precisely the same as that of reaction 3.
r1 )
k1XACVS KAAA KCYS KVAL ACV 1+ 1+ + + AAA CYS VAL KACV
(
)(
r2 )
k2XUPNSACV × O2 GSH ACV + Km 1 + Ki
(
(
r3 ) r4 )
)
))
k3XATIPN (IPN + Km,IPN)
(B.1)
(B.2)
(B.3)
k4XAT6APA × PoaCoa
(6APA × PoaCoa + PoaCoaKm,6APA-Poa + 6APAKm,Poa) (B.4)
r5 )
k5XATIPN × PoaCoa
(IPN × PoaCoa + Km,IPN-PoaPoaCoa + IPNKm,Poa) (B.5) r6 ) k66APA
(B.6)
r7 ) r3
(B.7)
r8 )
k8XATPenV (KPenV + PenV)
Material balance equations (B.14)-(B.21) describe the variation of the concentration of every species consid-1 31 ered inside the cell in mmol L-1 Cell h .
d(ACV) ) r1 - r2 - r9 - µACV dt
(B.14)
d(IPN) ) r2 - r3 - r5 - r10 - µIPN dt
(B.15)
d(6APA) ) (r3 - r4 - r6 + r8)Fcx dt
(B.16)
d(PenV) ) (r4 + r5 - r8)Fcx dt
(B.17)
d(ACVOUT) ) r9Fcx dt
(B.18)
d(IPNOUT) ) r10Fcx dt
(B.19)
d(8HPA) ) r6 dt
(B.20)
d(OPC) ) r7Fcx dt
(B.21)
In order to reproduce the results given by Pissara et al.,31 according to the experimental conditions in Jørgensen et al.,10,33 the following general fermentation equations were considered to describe a fed-batch system.
dV )F dt
(B.22)
d(xV) ) µxV dt
(B.23)
d(sV) µxV ) FsF dt Y
(B.24)
where µ is the specific growth rate given by the expression according to Contois kinetics proposed for penicillin fermentations:36
µ) (B.8)
r9 ) kACV(ACV - ACVOUT)
(B.9)
r10 ) kIPN(IPN - IPNOUT)
(B.10)
The intracellular specific activities of the enzymatic species are considered to vary during the fermentation,34 and that they are described by the following expressions:
d(XACVS) 1 )dt PACVS
(B.11)
d(XIPNS) 1 )dt PIPNS
(B.12)
d(XAT) 1 )dt PAT
(B.13)
µMs Ksx + s
(B.25)
The yield of the substrate is given as
1 1 m ) + Y Yo µ
(B.26)
and the maintenance coefficient m is given by
m)
mMs s + Ksm
(B.27)
The values for some of the parameters were taken from the literature, where these equations are extensively used for penicillin fermentation. Thus, Ks ) 0.006 g L-1 and Ksm ) 0.0001 g L-1.37,38 The rest were estimated to fit the specific data given by Jørgensen et al.10,33 so that µM ) 0.06 h-1, Yo ) 0.13 g g-1, and mM 0.02 g g-1
4708 Ind. Eng. Chem. Res., Vol. 37, No. 12, 1998
h-1. The parameters used in the equations presented in this appendix are shown in Table 4.31 Literature Cited (1) Cooney, C. L.; Acevedo, F. Theoretical conversion yields for penicillin synthesis. Biotechnol. Bioeng. 1977, 19, 1449. (2) Acevedo, F. Mass balancing: an effective tool for fermentation process optimization. CRC Crit. Rev. Biotechnol. 1987, 6, 309. (3) Voit, E. O. Optimization in integrated biochemical systems. Biotechnol. Bioeng. 1992, 40, 572. (4) Reagan, L.; Bogle, I. D. L.; Dunnill, P. Simulation and optimization of metabolic pathways. Comput. Chem. Eng. 1992, 16S, S237. (5) Reagan, L.; Bogle, I. D. L.; Dunnill, P. Simulation and optimization of metabolic pathways. Comput. Chem. Eng. 1993, 17, 627. (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) Delgado, J.; Liao, J. C. Inverse flux analysis for reduction of acetate excretion in Escherichia coli. Biotechnol. Prog. 1997, 13, 361. (9) Stephanopoulos, G.; Simpson, T. W. Flux amplification in complex metabolic networks. Chem. Eng. Sci. 1997, 52, 2607. (10) Jorgensen, H.; Nielsen, J.; Villadsen, J. Metabolic flux distributions in Penicilliun chrysogenum during fed-batch cultivations. Biotechnol. Bioeng. 1995a, 46, 117. (11) Schlosser, P; Bailey, J. An integrated modelling-experimental strategy for the analysis of methabolic pathways. Math. Biosci. 1990, 100, 87. (12) Sorribas, A.; Curto, R.; Cascante, M. Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis: Model definition and nomenclature. Math. Biosci. 1995a, 130, 25. (13) Sorribas, A.; Curto, R.; Cascante, M. Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis: Steady-state analysis. Math. Biosci. 1995b, 130, 51. (14) Sorribas, A.; Curto, R.; Cascante, M. Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis: Model validation and dynamic behavior. Math. Biosci. 1995c, 130, 71. (15) 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. (16) Savageau, M. Biochemical system analysis. i. some mathematical properties of the rate law for the component enzymatic reactions. J. Theor. Biol. 1969a, 25, 365. (17) 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. (18) Savageau, M. Biochemical system analysis. iii. dynamic solutions using a power-law approximation. J. Theor. Biol. 1970, 26, 215. (19) Savageau, M. Biochemical system theory: Operational differences among variant representations and their significance. J. Theor. Biol. 1991, 155, 509.
(20) Hatzimanikatis, V.; Bailey, J. E. Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log) linear kinetic models. Biotechnol. Bioeng. 1997, 54, 91. (21) Mauch, K.; Arnold, S.; Reuss, M. Dynamic sensitivity analysis for metabolic systems. Chem. Eng. Sci. 1997, 52, 2589. (22) Reder, C. Metabolic control theory: A structural approach. J. Theor. Biol. 1988, 135, 175. (23) Fell, D. Understanding of Control Metabolism; Portland Press: London, 1997. (24) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Ind. Eng. Chem. Res. 1994a, 33, 2111. (25) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a class of multistage dynamic optimization problems. 2. Problem with path constraints. Ind. Eng. Chem. Res. 1994b, 33, 2123. (26) Cuthrell, J. E.; Biegler, L. T. On the optimization of differential-algebraic process systems. AIChE J. 1987, 33, 1257. (27) Logsdon, J. S.; Biegler, L. T. Decomposition strategies for large-scale dynamic optimization problems. Chem. Eng. Sci. 1992, 47, 851. (28) Dunker, A. M. The decoupled direct method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys. 1984, 81, 2385. (29) Caracotsios, M.; Stewart, W. E. Sensitivity analysis of initial value problems with mixed odes and algebraic equations. Comput. Chem. Eng. 1985, 9, 359. (30) Schulz, A. R. Enzyme Kinetics: From Diastase to Multienzyme Systems; Cambridge University Press: Cambridge, U.K., 1994. (31) 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. (32) gPROMS User’s Guide, Version 0.1; Process Systems Enterprise Ltd.: London, 1996. (33) Jorgensen, 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. (34) 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. (35) Bryson, A. E.; Ho, Y.-C. Applied Optimal Control; Hemisphere Publishing Corp.: New York, 1975. (36) Bajpai, R. K.; Reuss, M. Evaluation of feeding strategies in carbon-regulated secondary metabolite production through mathematical modelling. Biotechnol. Bioeng. 1981, 23, 717. (37) Tayeb, Y. J.; Lim, H. C. Optimal glucose feed rates for fedbatch penicillin fermentation and computational results. Ann. N.Y. Acad. Sci. 1986, 469, 382. (38) San, K. Y.; Stepahnopoulos, G. Evaluation of feeding strategies in carbon-regulated secondary metabolite production through mathematical modelling optimization of fed-batch penicillin fermentation: A case of singular optimal control with state constraints. Biotechnol. Bioeng. 1989, 34, 72.
Received for review June 24, 1998 Revised manuscript received September 9, 1998 Accepted September 10, 1998 IE980410K
