Analysis of regulatory architectures in BST

PROGRESS IN NATURAL SCIENCEVol. 13 ,No.12 ,December 2003Analysis of regulatory architectures in BST*BIAN Fuping'2** . , LIANG Min' 2 and ZHAO Xueming'( 1. Department of M athematics , School of Science , Tianjin University ;2. Liu Hui Center for Applied M athematics , Nankai U niversity& Tianjin University ;3. School of Chemical Engineering , Tianjin University , Tianjin 300072 , China )Received April 11 , 2003 ; revised June 12 , 2003Abstractical systems is studied. After calculating the mixed-integer linear programming( MILP ) model given by Bailey et al. twice , the decisionmaking units( DMU ) and the prediction model of DEA are constructed , where the inputs are levels of manipulated parameters( enzyme )and outputs are concentrations of metabolites. W hen the metabolic net works are reconstructed , the data are obtained by calculating MILPframework twice and the optimal levels of the manipulated parameter at different regular loops are predicted , thus simplifying the calcula-tions of Bailey' :Keywords : metabolic reaction networks , MILP , the prediction model of DEA , eficient DMU.There are over one thousand kinds of enzymes inits optimal object.cells w hich catalyze various reactions and form a com-Introducing the constraints containing binaryplex reaction network. Due to the development ofbiochemistry and cellular physiology , researchers getvariables in S-system model we can get the mixed-in-to know the decomposition and synthesis pathways ofteger linear programming( MILP ) modef4]. Thesevarious components in cells and have a full under-new ly introduced constraints contain changes of a va-standing of controlling and regulation of these en-riety of enzy me regulatory architectures，reducing thezymes and their regulators. A large number of data ofamount of calculations. However ，w hen reactionenzyme' s dynamics have been accumulated throughpath ways become more complex , we need to solve thein vitro measurements. On this basis , through theMILP model many times. And the model becomesquantity analysis of the metabolic networks and de-very complicated as a result of introducing the newscription of the flux distribution of various metabolicconstraints. It is very difficult to solve the problempathways at different statuses in cells ，we can takeusing linear programming method , and the optimalsome improved measures and regulate their distribu-production rate cannot be predicted. In order to gettions and get more interesting products. Commonly ,the optimal net w ork architecture , this article predictsthe used methods include metabolic flux analysisthe production rate under different enzy me regular( MFA ), metabolic control analysis( MCA ) and bio-structures through utilizing the data envelopmentanalysis ( DEA ) and objective programming. It ischemical system theory( BST )1.based on the result of resolving the MILP twice andBST is an analytical method for metabolic net-the optimization of the metabolic net work .work developed in the 1970' s. On the basis of the1 Mixed-integer linear modeloptimization theory ，Voit et al.!23J found out onenet w ork architecture w hich optimizes the object func-We will consider that every reaction can be mod-tion. First of all , the relation between the reactionulated by any of the two metabolites ，X1 and X2，rates and their parameters ,e. g. concentrations of en-which will either inhihit or activate a reaction. Thiszyme , substrate , and reagents , is set up. Then thecons中国煤化工tulation of 12 regulato-variation range of the constraints is specified. The ob-:YHCNMHGjective function is the maximization of rate of produc-tion. By resolving the optimization problem ( S-sys-Four manipulated variables are considered : thetem model ) we can get the net work architecture withamount of the enzymes , P1 , P2 ,and P3 , that cat-# Supportgyt 书数ational Natural Science Foundation of China( Grant No. 20036010)** To whom cortp源dence should be ddressed. E-mail : fpbian @ public. tpt.tj. cnProgress in Natural ScienceVol.13 No.12 2003915strength Eij4=0.1 ,from metabolite j ;( ii) zi;j5s and Zij6 be equal to 1 , if reaction i isinhibited with strength Eε;js = - 0. 01 or activated .PPrwith strength Eij6=0.0l , from metabolite j.Then the S-system representation of the pathwayFig.1. A linear pathway with feedback inhibition. Bold solidis obtained4].lines denote reaction steps , dotted lines denote dependency on thecorresponding parameters Pl，dashed lines denote inhibition，andthin solid lines denote activation. .We can introduce a set of variables : .qi+qL=li(Pl)，L=1234,alyze the three reactions , and the amount of the effec-tor , P4 , that activates the first and the second reac-where qi denotes the logarithm of the reference valuetions of the pathway. Moreover ，for each loop we .of the parameter L , and qL denotes the logarithm ofwill consider Nreg = 6 alternative levels of regulatorythe factor by which the reference value is multipliedstrength and type of regulation :{-0.50.5 ,-0.1 ,to give the value pL. In the example studied here0.1 ,-0.01 0.01 }. We will allow only two regula-qi,=lnl=0tory loops to be active in the pathway .andqL=ln1=0Consider the manipulated parameter levels andthe regulatory structures that should be changed toMoreover，we introduce a set of binary vari-maximize the final product concentration X2 when theables，wI ,for which we will havefollowing conditions are satisfied :q' + wLqL = lr(Pl)， L= 1 234(i) The system is at a steady state ;and Z12 :vi= X;iB1口PP子,(ii) Xμ≤500;where v↑is net rate laws deseribing the processes(ii) V3≤10 ;andthat increase the concentration of metabolite 1，and( iv ) for the three enzy mes , only overexpressiong12 is kinetic orders.is considered- -that is ,Pl≥1( L=1 2 3 4) ,and upAt last , we introduce variables y; , tL , Sjm ,forto 10 times of their reference value.which we will haveWe will introduce the binary variables Zzjm andYj =Inxj，the parameters Ejm witht=WL9L'm=Sijm = EijmEijm"j ，respectively. These variables will be used in the de-j=l .. ,Nmet，scription of the steady-state equations after the loga-where Nreg is the number of the alternative strengthrithmic transformation. Thus we can write the MILPand types of regulation for each regulatory loop in themodet 41 :superstructure ,and Nrxn and N met are the numbers ofmaxy2the reactions and metabolites ，respectively , in thes.t.metabolic network. In this example , we have Nreg =6 , Nrn=3 and Nmet= 2. And for the binary vari-sjm +0.5y1 +! S2jim=1 j=1ables zjm and the parameters Ejm we let中国煤化工: In( 1/0.02),(i)zj1 and zqjz2 be equalto 1 ,if reaction i is in-JHCNMHG oo0.0y1 +"JS2jm -y2 -S3jm + 92hibited with strength εj;1 = - 0.5 or activated withstrength Eij2=0.5 ,from metabolite j ;+ t2 +2q4 +214- 9'3- t3 = lr( 2/0.02 ),y1≤Ir( 500),( ii)zij3 and Zzij4 be equal tol ,if reaction i is in-hibited w有万執据gth Eij3= -0.1 or activated withS12+ 9i+ 1+2q4+ 214≤Ir( 10),916Progress in Natural Science Vol.13 No. 122003(P≥1,2 DEA model and its efficiencyP2≥1，2.1 Fundamental definitionP3≥1,Ejmxy)-Sjm + mir( veim ryejm )xzimDefinition 1. The production possible set is{( x ry )loutput vector y can be obtained from input≥mir( y'ijm ry;eijm )EijmYj- Sijm + maxd( yjFijm ryjEijm )zjmi =1 23;x }.j=12;≤max( yjFijm ryjeijm )In this paper , when four enzyme expression lev-m=1p..6els are Pi ,P2 ,P3 ,P4 , respectively , the obtainableZijmmin( yEijm vyjeijm)- Sjm≤0intermediate product concentration is x1 and the finalSsjm - Zjmmax(( yjejm y;ijmn)≤0product concentration is x2.q- t1+ wqi≥qiDefinition2. If( x; y;) is an observed activity ，q1-t1+W1qi≤9ithen the referencesetis T={ x1 ry).. (xn ryn)}.l=1r.Awqi- I≤0Definition 3. When the relative increment per-tl-wIqi≤centage of input is more than that of the correspond-W1+W2+W3+W4≤1,ing output , the corresponding DMU of( x ,y )is de-creasing returns-to-scale ; W hen the relative increment之它之zjm = 2,percentage of input is less than that of correspondingm=l j=1 i=loutput，the corresponding DMU of( x ry )is increas-zjm≤l(i=123j=1 2).ing returns- to-scale ; w hen the relative increment percentage of input is equal to that of corresponding out-The best solution is found for this problem : .put, the corresponding DMU of( x y ) is constantx1= 500 ,x2= 125 ;returns-to-scale .PL= 1( L=1 3 4)，P2=2.236 ;01=1 i2.2 Fundamental DEA modelZ311= z321=1.Consider n DMU;( 1≤j≤n ). Their corre-Including the additional constraint in the modelsponding input vectors and output vectors areand solving the problem again , we find the secondx =(x1jr. xmj)">0,j=1r..n，best solution :yj =(yijr.. y)">0，j=1rnx1= 500，x2= 111.8 ;respectively. The intensity level of inputs and outputsPL=1( L=1 34)，P2= 10 ;are01= 10 ;0 =(01.. r0m)"Z311= 2Z212= 1.andu =( u1 r.. ,u,)respectively. The calculating of the intensity level ofAs indicated by the results above , after the con-inputs and outputs is based on certain laws.straint containing binary variable is introduced ，weinclude eight more linear constraints that w ill guaran-UyYkjtee the consistency between WLqL and t. ，ZzjmEjmnDjDefinition 4. h;=and Sijm，and thus the model becomes more compli-Dx;cated and the degree of difficulty in solving the modelis very high. Therefore , after getting the two resultsn ,is called the evaluation factor of efficiency of theabove，we do not introduce the MILP model again ,jth中国煤化工but analyze the new regulatory st ructure and predictthe optimal concentration based on the given data forMYHC N M H Gdicates that more out-the purpose of reducing the steps of iteration and con-puts are obtainable from fewer inputs.straint conditions .U sing Charnes-Cooper transformation , based onthe principle of duality of linear programming ，wecan construct a dual modetS] that has the non-Progress in Natural ScienceVol.13 No.12 2003917A rchimedes infinite small value ε :tion efficient for DEA.mir[θ-ε(e*s-+ eTs+)]= Vp。 ，Proof. See the footnote.n;:x;+ s~= 0x0，Theorem 1 gives a method for determining|j=whether DMU is effective when the inputs and out-( D。).t.{Ajj- st=yo,puts are regulated proportionally .λj≥0j=1rn,2.4 The prediction model of DEAs≥0,s*≥0.Consider how to predict the effective outputs ofThe purpose of model( De )is to get the maxi-the new DMU when there are a group of inputs andmum outputs with minimum inputs. When θ = 1 ,outputs of n DMU and an input of a new DMU. Itss~ =0 ,s+=0 ,( x0 ,yo) is efficient DMU ; whenalgorithm is as follows :θ<1 ( x0 ryo) is inefficient DMU of the DEA mod-Let Xmxn , Y;xn be the matrix composed of in-el. We can still produce the same outputs yo whileputs and outputs , respectively , Xo be the input of theconsuming fewer inputs.new DMU，Yo be the unknown outputs to be preThe limitation of this model is that it is not suit-dicted. First of all，we construct the following s pro-gramming problems :able for the production process that regulates the in-maxy;0i =12..sputs and outputs proportionally. To solve this prob-lem , Bian et al.6J have constructed a non-radial DEAλjij-yio=0,i=12.rs,model w hose inputs and outputs can be regulated pro-j=1portionally.s.t.\jxij ≤xi0，i =1 2r.,m2.3 N on- radial DEA model;=1,Now we introduce the non-radial DEA modelj= 1containing non-Archimedes infinite small value ε :λ;≥0，j =1 2r..n.1“1sminm二0;--β;-e(epSIp+ e's),Resolving these problems respectively ，we canget the ideal point of Yo :| 2x;xij+ Srp = 0:xi， i= 12r.. m，Y° =(yi0r20. ry).Secondly , we can set up a model containing weight :之入λjlij- s = Pryy，l=12..rs，max內Yyi0hereY; =一nyik2x=1，λj≥0，j=12..m,|j=1入jpij-yo=0,i=12mrs,Sp ,S≥0,0≤0;≤1，i=12..m,( D).t.xvqxij≤xi0，i=12rm,β≥1,l=12..s，=1,eIp=(1.1)°∈R”，le =(1.1)"∈R*.j≥0,j=12..n.The中国煤化工near programming mod-Theorem 1). The optimal solutions of this pro-el(Combining these twogramming problem are 0* ,i=1 .. ,m ,β° 1j=1 ,:JYHCNMHGgruup 1 outputs, which is the2.. 1s ,Sip ,S* ,and when 0:=1,i=1..,m ,output vector of the jth DMU.B;" =1,j=1.. ,s ,Sip=S* =O,DMU; is correc-1 ) Bian.The prediction model of DEA w ith undesirable outputs. Systems engineering-theory applications.918Progress in Natural Science Vol.13 No.12 20033 The DEA analysis in metabolic reactionmax y20networks[500λ1 + 500λ2- y10= 0,125λμ + 111.822- y20 = 0，For the metabolic reaction networks in Fig. 1 ,λ1+102≤1,we analyzed the optimal net w ork regulatory architec-( D2).t.2.236λ1+λ2≤1，tures by using the prediction model of DEA , consid-λ1+λ2=1，ering four enzy me expression levels P ,P2 ,P3，P4(λ1≥0入≥0as inputs and intermediate product concentration xand final product concentration x2 as outputs , andSolving the model( D1 ) and model( D2 )( calcu-constructed DMU by using prediction model of DEA.late them by using the linprog function in Matlab ) ,We determined if the reaction pathway in Fig.1 is op-wegettheoptimalpointofXo:Xg=(500,125).timal according to the efficiency of DMU and thenconstructed4 x 2 input matrix and 2 X 2 output ma-The second step is constructing the model con-trix :taining weight :「P11 P12_Y10_Y20 .P21 P22[Xu X12]maxP4x2=X2x2=( 500+500) 一( 125+ 111.8)P31 P32X21 X22J[500λ1 + 500λ2- y1o=0 ,125λ1+ 111.8λ2- y2o=0 ,The elements of the ith column in P4x2 are thelevel of enzyme of the ith reaction pathway ,i=1 2.(D3)3.t.λ1+ 10λ2≤1 ，| 2.236λ1 + λ2≤1 ,The first row and the second row elements of the ithλ1+ λ2=1 ,column in X2x2 are the intermediate product concen-lλi≥0小2≥0.tration and final product concentration of the ith re-action pathway , respectively .The optimal solution of the model( D3 )is( 500 ,Substituting the data of Pr，P2 , P3 , P4 and125 ).X1，X2 into the matrixes above after solving theIn summary ,it is clear that the optimal produc-MILP model twice , we havetion of the third reaction pathway is( 500 ,125 ). This「1107conclusion is the same as the optimal production with2.236:[ 500 500P4x2=| 11 | ' A2x2=l125 111. 8Jmultiple construction MILP in Ref. [4 ]. The pro-duction calculated from S-system model is , however,L1( 100 ,5 ). Obviously , this pathway is not optimal.Now , we predict the optimization of the thirdThe prediction DEA model calculation used in this pa-reaction pathway by constructing DM Uo , letting theper is simpler than MILP. For the reaction in w hichlevel of 4 kinds of enzymes( P: ,P2 ,P3 ,P4)=(1 ,1，more than three pathways exist , we can reset the en-1 5 2 ) be the inputs ,and predicting its optimal out-zyme expression levels and regulatory architectureslist the prediction DEA model several times , and preputs.dict the new reaction pathway ，and iterate like thisAccording to the steps of the predicted algo-until we obtain an optimal reaction pathway .rithm , the first step is setting up s programming4Conclusionproblems as follows :maxy1o心he eficient DMU by us-( 500λ1 + 500λ2- y1o= 0 ,ing中国煤化工me )and predicted out-125λ1 + 111.8λ2- y20= 0，puts,RYHCNM H Gts). It resolves the pre-λ+10λ2≤1，dicted optimal outputs and compares them with the .(DI ).t.2.236λ1 + λ2≤1 ,actual value of the outputs. Through observing the .λ1+λ2=1,optimization of the reaction , we confirmed w hetheror not to regulate manipulated parameters and recon-1≥022≥0，structed predicted DEA model. Due to the simplicityProgress in Natural ScienceVol.13 No.12 2003919of predicting DEA model and the possibility of calcu-3 Voit, E. 0. Optimization in integrated biochemical systems.lating them using the given procedure , the amount ofBiotechnol. Bioeng. , 1992 ,40(5 ):572.Vassily, H. et al. Optimization of regulatory architectures incalculation is reduced dramatically , indicating thatmetabolic reaction net works. Biotechnology and Bioengineeringthis method is feasible .1996 ,52 :485.i Sheng,Z. H. et al. Theory , Method and Applying of DEA( inReferencesChinese ). Beijing : Science Press , 1996 :30 ~ 36.6 Bian ,F. P. et al. The fficiency valuation of DM Us with1 Zhao ,X. M. et al. Metabolic engineering : its history ,present andable outputs. Transaction of Tianjin University ( in Chinese ),future. Chemical Engineering( in Chinese ), 1997 :1319.002 ,84):295 ~ 298.2 Torres, N. V. et al. Optimization of nonlinear biotechnologicalprocesses with linear programming : application to citric acid pro-duction by Asperilus niger. Biotechnol. Bioeng. , 1996 ,49 3):中国煤化工MYHCNMHG.