Modelling and pathway identification involving the transport mechanism of a complex metabolic system in batch culture

Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Modelling and pathway identification involving the transport mechanism of a complex metabolic system in batch culture Jinlong Yuan a,b,⇑, Xu Zhang a, Xi Zhu a, Enmin Feng a, Hongchao Yin b, Zhilong Xiu c a

School of Mathematical Science, Dalian University of Technology, Dalian, Liaoning 116024, PR China School of Energy and Engineering, Dalian University of Technology, Dalian, Liaoning 116024, PR China c School of Environmental and Biological Science and Technology, Dalian University of Technology, Dalian, Liaoning 116012, PR China b

a r t i c l e

i n f o

Article history: Received 2 June 2013 Received in revised form 26 September 2013 Accepted 20 October 2013 Available online 1 November 2013 Keywords: Pathway identification Complex metabolic system Biological robustness Constraint transcription Parallel MPSO

a b s t r a c t The bio-dissimilation of glycerol to 1,3-propanediol (1,3-PD) by Klebsiella pneumoniae (K. pneumoniae) can be characterized by a complex metabolic system of interactions among biochemical fluxes, metabolic compounds, key enzymes and genetic regulation. In this paper, in consideration of the fact that the transport ways of 1,3-PD and glycerol with different weights across cell membrane are still unclear in batch culture, we consider 121 possible metabolic pathways and establish a novel mathematical model which is represented by a complex metabolic system. Taking into account the difficulty in accurately measuring the concentration of intracellular substances and the absence of equilibrium point for the metabolic system of batch culture, the novel approach used here is to define quantitatively biological robustness of the intracellular substance concentrations for the overall process of batch culture. To determine the most possible metabolic pathway, we take the defined biological robustness as cost function and establish an identification model, in which 1452 system parameters and 484 pathway parameters are involved. Simultaneously, the identification model is subject to the metabolic system, continuous state constraints and parameter constraints. As such, solving the identification model by a serial program is a very complicated task. We propose a parallel migration particle swarm optimization algorithm (MPSO) capable of solving the identification model in conjunction with the constraint transcription and smoothing approximation techniques. Numerical results show that the most possible metabolic pathway and the corresponding metabolic system can reasonably describe the process of batch culture. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction The microbial conversion of glycerol by Klebsiella pneumoniae (K. pneumoniae) to 1,3-propanediol (1,3-PD) in batch culture is of interest to industry because of its environmentally safe, high region specificity, cheaply available feedstock, and relatively high theoretical molar yield [1–5]. Therefore, since the 1980s, many a computational model has been established to describe the process of batch culture. Based on the proposed and modified model by Zeng et al. [6,7], researchers investigated optimal control [8], multistage model [9], parameter identification [10] and stochastic model [11] in batch culture.

⇑ Corresponding author at: School of Mathematical Science, Dalian University of Technology, Dalian, Liaoning 116024, PR China. Tel.: +86 41184708351x8025; fax: +86 41184708354. E-mail address: [email protected] (J. Yuan). 1007-5704/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cnsns.2013.10.021

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

2089

The concentrations of intracellular substances are almost always ignored in any of early batch culture references mentioned above. This is obviously a serious limitation, as intracellular substances arising in real-world applications really exist. In the context of less information about intracellular experiment data, the quantitative description of biological robustness becomes a feasible method to evaluate the validity of the computational concentrations of intracellular substances. Robustness is one of the fundamental characteristics of biological system and it allows a system to sustain its functions in spite of internal and external perturbation [12–14]. Marhl and Perc [15] defined the robustness of the system via a given parameter and a dependant variable, like for example in their case peak and plateau value. Usually, the robustness of the system is evaluated by the sensitivity analysis [16] and jitter [17]. Perc et al. [18,19] studied only the frequency robustness of the system since the amplitude of the new spike is the same as the amplitude of the original spikes. This point of view, which has been observed for a wide variety of experiments [20,21], is being gradually accepted by experts in the field of systems biology. The biological robustness has been quantitatively defined in the approximately stable state of continue culture [22,23]. None of the above literatures, however, quantitatively define the biological robustness of batch culture owing to lack of equilibrium point. Up to now, the biological robustness defined quantitatively is seldom found in batch culture. In 2008, Sun et al. [24] firstly proposed a novel mathematical model to describe the concentration changes of extracellular and intracellular substances. This model assumed that glycerol passes cell membrane by passive diffusion coupled with facilitated transport and 1,3-PD is transported by passive diffusion. To the best of our knowledge, since there is a small number of reports dealing with the transport ways of glycerol and 1,3-PD of batch culture and it is still not exactly known, the reliability of the model cannot be guaranteed [25]. In this paper, taking into account three possible transport ways of glycerol and 1,3-PD with different weights across cell membranes (passive diffusion, active transport or passive diffusion coupled with active transport), we give 121 possible metabolic pathways and establish a novel mathematical model, which is represented by a complex metabolic system. Some properties of the metabolic system are discussed. Similar to literature [18,19], inspired by the qualitative description of the biological robustness given by Kitano [12–14] and h-stability [26], we first put forward a quantitative definition of biological robustness of intracellular substance concentrations for the batch culture because of less information concerning intracellular substances. The proposed biological robustness, which is defined for the overall process of batch culture due to the absence of equilibrium point, is different from quantitative description of biological robustness of continue culture in the approximately stable state [22,23]. Taking the proposed biological robustness as a cost function, we establish an identification model involving the metabolic system together with subject to continuous state constraints and parameter constraints. The constraint transcription and smoothing approximation techniques [27] are applied to deal with the continuous state constraints in the identification model. The case that the model includes 1452 system parameters and 484 pathway parameters makes solving the identification model for a serial program a very complicated task. With this in mind, we construct a parallel migration particle swarm optimization algorithm to determine the most possible metabolic pathway and corresponding optimal parameters under various initial conditions. Applying the proposed algorithm, in which 9979188 numerical computations of differential equations are involved, to the identification model on 1800 PC-cluster Server, we obtain that both glycerol and 1,3-PD pass cell membranes by active transport and passive diffusion with different weights. The remainder of this paper is organized as follows. In Section 2, a complex metabolic system is formulated and some important properties are proved. In Section 3, the identification model is proposed via biological robustness. In Section 4, a parallel particle swarm optimization algorithm based on particle migration (MPSO) is constructed to solve the identification model. In Section 5, numerical results are presented. In Section 6, we draw the conclusions and trace the direction for future works.

2. Complex metabolic system The fermentations of glycerol cover both extracellular and intracellular environments, which are linked by the transports of substrate (glycerol) and products (1,3-PD) across cell membranes. As shown in Fig. 1 [24], the transport ways of glycerol and 1,3-PD across the membrane have not been observed in experiments yet. Abbreviations: A, active transport; P, passive diffusion; AP, passive diffusion and active transport; GTW, transport ways of glycerol; PDTW, transport ways of 1,3-PD. Taking three possible transport ways of glycerol and 1,3-PD across cell membranes{A,P,AP} into consideration, let D :¼ fð1; sÞ 2 ½0; 12 g be the set of possible metabolic pathways, where 1ands are pathway parameters. Since 1 and s are continuous variables, we can obtain infinite metabolic pathways. The transport ways of glycerol corresponding to parameter 1 are as follow:

8 > < 1 :¼ 0; if GTW is P; 1 :¼ 1; if GTW is A; > : 1 2 ð0; 1Þ; if GTW is AP: The transport ways of 1,3-PD corresponding to parameter

s are as follow:

2090

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

Fig. 1. Anaerobic metabolism pathways of glycerol fermentation [24]. Abbreviations:GDHt, glycerol dehydratase; 3-HPA, 3-hydroxypropionaldehyde; 1,3PD, 1,3-propanediol; PDOR, 1,3-PD oxydoreductase; GDH, glycerol dehydrogenase; DHA, dihydroxyace-tone; DHAK, dihydroxyacetone kinase; DHAP, dihydroxyacetonephosphate; HAc, acetic acid; ETOH, ethnol; TCA, TCA cycle.

8 > < s :¼ 0; if PDTW is P; s :¼ 1; if PDTW is A; > : s 2 ð0; 1Þ; if PDTW is AP: It follows from the above description that the main goal of this paper is to find the possible metabolic pathway. In the light of the factual experiment, the following conditions will be in force throughout the rest of this paper, whether explicitly mentioned or not. H1. No medium is pumped inside and outside the bioreactor in the process of batch fermentation. H2. The concentrations of reactants are uniform in reactor, while time delay and nonuniform space distribution are ignored. Let xðtÞ ¼ ½x1 ðtÞ; . . . ; x8 ðtÞT 2 R8 be the continuous state vector, where x1 ðtÞ; . . . ; x8 ðtÞ denote the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol, intracellular glycerol, 3-HPA and intracellular 1,3-PD at time t 2 ½0; tf , respectively. ½0; t f   Rþ is the interval of reaction time, and x0 2 R8 , is the initial state. To simplify notation, we denote u1;s as u; xi ðtÞ as xi ; i 2 I8 , if there is no confusion. Under the assumptions (H1) and (H2), for given 1; s 2 ½0; 1, mass balance of biomass, substrate and products in batch culture can be described as the following complex metabolic system in literature [10], denoted by CMNð1; sÞ



_ xðtÞ ¼ hðx; u1;s ; 1; sÞ; t 2 ½0; t f  xð0Þ ¼ x0 :

ð1Þ

12 where u1;s ¼ ðu11;s ; u21;s ; . . . ; u12 1;s Þ 2 R denotes the system parameters to be estimated. The right hand side of system (1) is of T the form hðx; u1;s ; 1; sÞ :¼ ðh1 ðx; u1;s ; 1; sÞ; . . . ; h8 ðx; u1;s ; 1; sÞÞ with the components [10,24] defined as

8 h1 ðx; u; 1; sÞ ¼ lx1 ; > > > > > h2 ðx; u; 1; sÞ ¼ q2 ð1; sÞx1 ; > > > > > hj ðx; u; 1; sÞ ¼ qj ð1; sÞx1 ; j ¼ f3g > > > > < hj ðx; u; 1; sÞ ¼ qj x1 ; j ¼ f4; 5g h i 2 > h6 ðx; u; 1; sÞ ¼ u19 1 x46:1633x þ ð1  1Þu10 ðx2  x6 ÞNRþ ðx2  x6 Þ  q20  lx6 ; > > þ2:7467 2 > > > > x 45:9992x > > h7 ðx; u; 1; sÞ ¼ u11 0:53ð1þ 6x7 Þþx  0:14þx ð1þ 7x7 Þ  lx7 ; > 7 6 > 185:242 1:3341 > > > 6:9401x8 45:9992x7 : h8 ðx; u; 1; sÞ ¼ 0:14þx  s  ð1  sÞu12 ðx8  x3 ÞNRþ ðx8  x3 Þ  lx8 : x7 x8 þ26:6321 ð1þ Þ 7

1:3341

The specific substrate consumption rate and the specific product formation rate can be expressed as follow:

ð2Þ

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

8 l 2 > q20 ¼ 3:1453 þ 0:0066 þ x36:8436x ; > 2 þ16:2526 > > Q > x 5 0:67x2 j > > > > l ¼ x2 þ0:28 j¼2 ð1  xj Þ; > > < x2 q2 ð1; sÞ ¼ 1u1 x2 þ3:87 þ ð1  1Þu2 ðx2  x6 ÞNRþ ðx2  x6 Þ; > 5 > 8 > q3 ð1; sÞ ¼ su3 x xþu 4 þ ð1  sÞu ðx8  x3 ÞN Rþ ðx8  x3 Þ; > 8 > > > 6 > > q4 ¼ 1:0522 þ lu ; > > : q5 ¼ u7 þ lu8 :

2091

ð3Þ

Different transport ways of glycerol and 1,3-PD correspond to the different metabolic pathways, and the different operation condition in actual experiment correspond to the different parameters [28]. Accordingly, the expressions of q2 ð1; sÞ; q3 ð1; sÞ; h6 ðx; u; 1; sÞ; h8 ðx; u; 1; sÞ are different as the possible transport ways are active transport, passive diffusion or passive diffusion coupled with active transport with 1; s 2 ½0; 1. In addition, the indicator function

 NRþ ðhÞ ¼

1; h > 0; 0;

h 6 0;

means that the structures of equations which contain this term are variant. As the same with the previous work [10], u0 :¼ ½51:4459; 101:122; 144:066; 3:7595; 85:2486; 53:3623; 2:8003; 12:8892; 8:2813; 4543:09; 38:7706; 71:8511T 2 R12 . The admissible set of system parameters u 2 R12 is defined as:

U ad :¼ ½u ; u  ¼

12 Y

½uj ; uj  ¼

j¼1

12 Y ½u0j  ju0j j=2; u0j þ ju0j j=2  R12 : j¼1

In the light of literature [24], x ¼ ½0:0001; 0:01; 0; 0; 0; 0; 0; 0T ; x ¼ ½15; 2039; 939:5; 1026; 360:9; 2039; 300; 939:5T . Under anaerobic conditions at 37 °C and pH¼ 7:0, since each component of the state variable xðtÞ represents a certain substance concentration, the concentrations of biomass, glycerol and products should be restricted in a certain range according to the practical production. Given xj ; j 2 I8 , is the critical concentration of xj ðtÞ (the maximum of xj ðtÞ) and xj P 0; j 2 I8 , is the lower bound of xðtÞ. So the admissible range W 0 of xj ðtÞ is:

xðtÞ 2 W 0 :¼ ½x ; x  ¼

8 Y

½xj ; xj   R8þ :

ð4Þ

j¼1

3. Properties of the complex metabolic system The purpose of this section is to demonstrate some properties of the complex metabolic system. From (2) and (3), we have the following Property 1: Property 1. For given 1;

s 2 ½0; 1,

(I) 8x 2 W 0 ; hðx; u; 1; sÞ is continuous in u on U ad . (II) 8u 2 U ad ; hðx; u; 1; sÞ is Lipschitz continuous in x on W 0 . (III) 8u 2 U ad ; hðx; u; 1; sÞ satisfies the linear growth condition, namely, 8x 2 W 0 ; there exist f; # > 0, such that khðx; u; 1; sÞk 6 fkxk þ #, where k  k is a Euclidean norm.

Proof (I) The fact that every component function of hðx; u; 1; sÞ is continuous in u on U ad from the fact xðtÞ is continuous in t on ½0; t f . Let us note that hðx; u; 1; sÞ is vector valued function, hence, hðx; u; 1; sÞ is continuous in u on U ad . (II) 8x; y 2 W 0 ; hj ðx; u; 1; sÞ; j 2 I8 , is continuously differentiable in x on W 0 . It follows from Mean Value Theorem that

hj ðx; u; 1; sÞ  hj ðy; u; 1; sÞ ¼ h5x hj ðhx þ ð1  hÞy; u; 1; sÞ; ðx  yÞi; j 2 I8 ; where h 2 ½0; 1. It follows from the Cauchy–Schwartz inequality that

jhj ðx; u; 1; s; kÞ  hj ðy; u; 1; sÞj ¼ jh5x hj ðhx þ ð1  hÞy; u; 1; sÞ; ðx  yÞij 6 k5x hj ðhx þ ð1  hÞy; u; 1; sÞk  kx  yk: Since W 0  R8 is a compact set and hj ðx; u; 1; sÞ; j 2 I8 , is continuously differentiable in x on W 0 , we can deduce

k5x hj ðhx þ ð1  hÞy; u; 1; sÞk 6 Lm :¼ max j5x hj ðhx þ ð1  hÞy; u; 1; sÞj; x;y2W 0

2092

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

which shows that jhj ðx; u; 1; sÞ  hj ðy; u; 1; sÞj 6 Lm kx  yk. Hence, hj ðx; u; 1; sÞ; j 2 I8 , is Lipschitz continuous in x on W 0 . Combining the above results and hðx; u; 1; sÞ is a vector valued function, the desired result follows. (III) Now, it follows from (3) that

8 l  0:67; > > > > l > > j þ 36:8436; q20  3:1453 þ j 0:0066 > > > > < q ð1; sÞ 6 j1j  ju1 j þ jð1 þ 1Þj  ju2 j  jðx2 þ x6 Þj; 2 3 5 > q > 3 ð1; sÞ 6 jsj  ju j þ jð1 þ sÞj  ju j  jðx8 þ x3 Þj; > > > > 6 > q4  1:0522 þ jlj  ju j; > > > : q5 6 ju7 j þ jlj  ju8 j; which shows that

8 jh1 ðx; u; 1; sÞj 6 jlj  jx1 j; > > > > > jhj ðx; u; 1; sÞj 6 jqj ð1; sÞj  jx1 j; j 2 f2; 3g > > > > > < jhj ðx; u; 1; sÞj 6 jqj j  jx1 j; j 2 f4; 5g > jh6 ðx; u; 1; sÞj 6 j u19 j  ð46:1633  j1j þ jð1 þ 1Þj  ju10 j  ðjx2 j þ jx6 jÞ þ jq20 jÞ þ jlj  jx6 j; > > > > > > jh7 ðx; u; 1; sÞj 6 ju11 j þ 45:9992 þ jlj  jx7 j; > > > : jh8 ðx; u; 1; sÞj  45:9992 þ 6:9401  jsj þ jð1 þ sÞj  ju12 j  ðjx8 j þ jx3 jÞ þ jlj  jx8 j:

ð5Þ

The fact that u is bounded on U ad follows from the fact that U ad is bounded and that u 2 U ad . It can be observed that all of the coefficients of jxj jðj 2 I8 Þ and the remaining terms are constants in the right-side of inequalities (5). Let K 1 be apmaximal ffiffiffi value p offfiffiffithe coefficients of jxj jðj 2 I8 Þ, and K 2 be the maximum of the remaining terms. At the last, letting f :¼ 2 2K 1 and # :¼ 2 2K 2 , we obtain the desired result. h Armed with the classical theory of differential equation and Property 1 [29–31], we have the following property: Property 2. For given 1; s 2 ½0; 1; 8x0 2 W 0 ; 8u 2 U ad , the solution xðtÞ of CMNð1; sÞ, denoting by xðt; x0 ; u; 1; sÞ, is existence and uniqueness; and xðt; x0 ; u; 1; sÞ is continuous with respect to the parameter vector u 2 U ad and the initial value x0 2 W 0 .

4. Biological robustness index The purpose of this section is to establish a reliable criterion to determine the most possible metabolic pathway. Let IN :¼ f1; 2; . . . ; Ng be the serial number set of experiments in batch culture, where N is the total experiment times. xl0 2 W 0 is the initial state of the l 2 IN experiment. For given 1; s 2 ½0; 1; 8l 2 IN , let xð; xl0 ; u; 1; sÞ be the solution of CMNð1; sÞ and xj ð; xl0 ; u; 1; sÞ, be the jth 2 I8 component of xð; xl0 ; u; 1; sÞ. For given 1; s 2 ½0; 1, we define the following set:

S0 ð1; sÞ :¼ fxð;xl0 ;u; 1; sÞjxð;xl0 ;u; 1; sÞ 2 W 0 ; is the solution of CMN ð1; sÞ with initial state xl0 for u 2 U ad ; 8l 2 IN g:

ð6Þ

The set of feasible parameter vectors is defined as

U S0 ð1; sÞ :¼ fu 2 U ad jxðt; xl0 ; u; 1; sÞ 2 S0 ð1; sÞ; 8t 2 ½0; tf ; 8l 2 IN g:

ð7Þ

In view of the factual experiments, we can assume that. H3. For given 1;

s 2 ½0; 1, the sets S0 ð1; sÞ and U S0 ð1; sÞ are nonempty, respectively.

It follows from the compactness of U ad and W 0 , together with the continuity of the solutions with respect to u 2 U S0 ð1; sÞ that we will have frequent occasion to obtain the following result: Property 3. Under the assumption H3, given 1; s 2 ½0; 1, the sets S0 ð1; sÞ and U S0 ð1; sÞ are compact, respectively. 4.1. Relative error In order to evaluate the reliability of system (1), we should give an evaluation criterion for this system. We can only obtain experiment data of extracellular substances. In a general way, the computational results are respected to be in

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

2093

accordance with experimental data. Since the feeding of alkali whose purpose is to maintain its pH value at 7.0 or so leads to the inaccuracy of the extracellular concentrations of acetic acid and ethanol, this paper is concerned with the relative error between experimental data and computational values of the first three substances. Let yl ðtÞ ¼ ðyl1 ðtÞ; . . . ; yl5 ðtÞÞ 2 C 1 ð½0; t f ; R5 Þ; l 2 IN , be the vector function fitted by experiment data with initial state xl0 2 W 0 , namely, it’s components denote the concentration of every extracellular substances at time point t 2 ½0; tf . Definition 1. Given 1; s 2 ½0; 1, the relative error between computational concentrations of extracellular components xj ðt; xl0 ; u; 1; sÞ; j 2 I3 and experimental concentrations ylj ðtÞ; j 2 I3 can be defined as

SSDðu; 1; sÞ :¼

N X 3 Z tf ½x ðt; xl ; u; 1; sÞ  yl ðtÞ2 X 1 j 0 j dt: 2 3  N  tf l¼1 j¼1 0 ½yl ðtÞ

ð8Þ

j

In this paper, the relative error of the extracellular substance concentrations is limited in ½0; n; n is the tolerance of relative error SSDðu; 1; sÞ. By (6)–(8), given 1; s 2 ½0; 1, we define

U w ð1; sÞ :¼ fuju 2 U S0 ð1; sÞ; SSDðu; 1; sÞ 6 ng:

ð9Þ

If the set U w ð1; sÞ is empty, we argue that CMNð1; sÞ is very poor in characterizing the actual fermentation process and let its optimal biological robustness be þ1. Otherwise, we will calculate biological robustness of CMNð1; sÞ. In view of the factual fermentation, we make the following assumption. H4 . 9 1; s 2 ½0; 1, such that the set U w ð1; sÞ is nonempty. Due to the compactness of U S0 ð1; sÞ by Property 3, we will have frequent occasion to obtain the following result: Property 4. Under the assumption H4, for some 1; s 2 ½0; 1, the sets U w ð1; sÞ is compact. 4.2. Biological robustness The analytical solution to system (1) is not easy to find, therefore numerical approach is introduced to overcome the difficulty. When considering intracellular substances, how to judge the reliability of numerical solution to system (1) is the main purpose of the section because of lack of intracellular data. Taking into account the basic features of biological systems, which might be helpful to address the issue. One of the fundamental and ubiquitously observed characteristics of biological systems is robustness. It is a property that allows a system to maintain its functions in spite of external and internal perturbations [12–14]. The biological robustness has been quantitatively defined in the approximately stable state of continue culture [22,23]. Marhl and Perc [15] defined the robustness of the system via a given parameter and a dependant variable, like for example in their case peak and plateau value. Usually, the robustness of the system is evaluated by the sensitivity analysis [16] and jitter [17]. The above methods are improper to define biological robustness of batch culture because the system of characterizing the process of batch culture is the absence of equilibrium point. Within such system of characterizing the process of batch culture, the word ‘robust’ is frequently used, rarely defined, and never quantified. Hence, it is difficult to compare the proposed definition with other definitions of biological robustness in batch culture. The definition of biological robustness in batch culture should meet the case – the smaller the changes in state variables provoked by parameters small shift for the overall process of batch culture, the more robust the system is. Perc et al. [18,19] studied only the frequency robustness of the system since the amplitude of the new spike is the same as the amplitude of the original spikes. With this in mind, inspired by the qualitative description of the biological robustness and h-Stability [26], similar to literature [18,19], we give the quantitative definition of biological robustness, which is completely different from quantitative description of biological robustness of continue culture in the approximately stable state, for the overall process of batch culture to describe the variation about the parameter disturbance u0  u due to the continuity of xð; xl0 ; u; 1; sÞ; l 2 IN with respect to u 2 U w ð1; sÞ as follow. Definition 2. Given 1; is defined as:

s 2 ½0; 1, for u 2 U w ð1; sÞ, the biological robustness of system parameter of the solution of CMNð1; sÞ

N Z 1 X Jðu; 1; sÞ :¼ /ðu0 Þ N  t f l¼1 U

R tf 0

kxðt;xl0 ;u0 ;1;sÞxðt;xl0 ;u;1;sÞk kxðt;xl0 ;u;1;sÞk

ku0

 uk

dt dðu0 Þ;

ð10Þ

where U  U w ð1; sÞ is the perturbation space with u, and /ðu0 Þ is taken as probability density function of uniform distribution in u0 on U. Given 1; s 2 ½0; 1, we can conclude that the smaller the value of Eq. (10) is, the more robust of the system (1) is, i.e., given u1 2 U w ð1; sÞ; u2 2 U w ð1; sÞ, the system (1) with regard to its parameter vector u1 against a set of perturbations in Bðu1 ; dÞ can be said to be more robust than with regard to its other parameter vector u2 against a set of perturbations in Bðu2 ; dÞ, when

2094

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

Jðu1 ; 1; sÞ < Jðu2 ; 1; sÞ: 5. Pathway identification and parallel MPSO algorithm In this section, an identification model together with its approximate model for the complex metabolic system is established and a parallel particle swarm optimization algorithm based on particle migration (MPSO) is constructed to solve it. 12 Given 1; s 2 ½0; 1, let u1;s ¼ ðu11;s ; u21;s ; . . . ; u12 1;s Þ 2 R be the identified parameter vector of CMNð1; sÞ. 5.1. Pathway identification model The pathway identification model of CMNð1; sÞ; 1; s 2 ½0; 1, can be viewed as:

ModelðOIPÞ : min :Jðu1;s ; 1; sÞ s:t: u1;s 2 U w ð1; sÞ

1; s 2 ½0; 1 It follows from the above description that the following theorem is to show the identifiability of the model OIP. Theorem 1. For the model OIP, given 1; s 2 ½0; 1, there exists u1;s 2 U w ð1; sÞ such that

Jðu1;s ; 1; sÞ 6 Jðu1;s ; 1; sÞ; 8u1;s 2 U w ð1; sÞ:

ð11Þ

Proof. Given 1; s 2 ½0; 1, it follows from Properties 2–4 that Jðu1;s ; 1; sÞ is a continuous function on the compact set U w ð1; sÞ. Hence, the theorem must be true. h In above model OIP, because of the coexistence of two types of parameters, i.e., system parameters u1;s 2 U w ð1; sÞ, and pathway parameters 1; s 2 ½0; 1, it is difficult to solve it directly. In addition, system parameters and pathway parameters are independent, so we can split the OIP into two subproblems. Firstly, given 1; s 2 ½0; 1, the identification subproblem corresponding to CMNð1; sÞ, denoted by OIP1 ð1; sÞ, can be formulated as

ModelðOIP1 ð1; sÞÞ : J 1 ðu1;s ; 1; sÞ :¼ min :Jðu1;s ; 1; sÞ s:t: u1;s 2 U w ð1; sÞ Since 1; s 2 ½0; 1, are given, there is only continuous parameter vector u1;s 2 U w ð1; sÞ in OIP1 ð1; sÞ. Secondly, the identification subproblem concerning 1; s 2 ½0; 1, denoted by OIP2 , can be formulated as

ModelðOIP2 Þ : J 2 ðu1 ;s ; 1 ; s Þ :¼ min J 1 ðu1;s ; 1; sÞ; 1;s2½0;1









ðu1 ;s ; 1 ; s Þ :¼ arg min J 1 ðu1;s ; 1; sÞ: 1;s2½0;1

A mathematical programming problem with a finite number of variables and infinitely many constraints is called a semiinfinite programming problem [40]. Clearly, in view of the definition of semi-infinite programming problem and the above description, model OIP, OIP1 ð1; sÞ and ðOIP2 Þ can be viewed as semi-infinite programming problems because of the continuous state inequality constraint (4) and 1; s 2 ½0; 1. 5.2. Approximate model When we intend to solve the problem OIP, our main difficulty lies in solving the subproblem OIP1 ð1; sÞ; 1; s 2 ½0; 1. In nature, given 1; s 2 ½0; 1; OIP1 ð1; sÞ is a constrained optimization problem. In addition, based on the above description, we see that OIP1 ð1; sÞ is a semi-infinite programming problem. An efficient algorithm, in which the so-called constraint transform and local smoothing techniques are involved, can address the issue [27]. Given l 2 IN ; 1; s 2 ½0; 1, for the model OIP1 ð1; sÞ, it is difficult to directly judge whether the constraint condition u1;s 2 U w ð1; sÞ holds or not. In practice, the essential difficulty lies in the judgement of the constraint

xðt; xl0 ; u1;s ; 1; sÞ 2 S0 ð1; sÞ; 8t 2 ½0; t f :

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

2095

To surmount these difficulties, let

g j ðt; l; u1;s ; 1; sÞ :¼ xj  xj ðt; xl0 ; u1;s ; 1; sÞ; g 8þj ðt; l; u1;s ; 1; sÞ :¼ xj ðt; xl0 ; u1;s ; 1; sÞ  xj ; j 2 I8 : The constraint xðt; xl0 ; u1;s ; 1; sÞ 2 S0 ð1; sÞ is equivalently transcribed into

Gðl; u1;s ; 1; sÞ ¼ 0;

ð12Þ

P16 R tf

where Gðl; u1;s ; 1; sÞ :¼ j¼1 0 minf0; g j ðt; l; u1;s ; 1; sÞgdt. However, Gðl; u1;s ; 1; sÞ is non-smooth in u1;s 2 U w ð1; sÞ. Consequently, standard optimization routines would have difficulties with this type of equality constrains. we replace (12) with 16 Z X

e e;c ðl; u1;s ; 1; sÞ :¼ c þ G

j¼1

where

tf

0

ue;j ðt; l; u1;s ; 1; sÞdt P 0;

ð13Þ

e > 0; c > 0 and 8 g ðt; l; u1;s ; 1; sÞ; if > > < j 2 ue;j ðt; l; u1;s ; 1; sÞ :¼  ðgj ðt;l;u1;se;1;sÞeÞ ; if > > : 0; if

g j ðt; l; u1;s ; 1; sÞ < e;  e 6 g j ðt; l; u1;s ; 1; sÞ 6 e; g j ðt; l; u1;s ; 1; sÞ > e:

Therefore, OIP1 ð1; sÞ; OIP2 and OIP can be approximated by a sequence of nonlinear programming models fOIP1;e;c ð1; sÞg; fOIP2;e;c g and fOIPe;c g defined by replacing constraint (12) with (13), respectively. As is shown in [27], the following theorem shows that the solution of the corresponding problem fOIP1;e;c ð1; sÞg; fOIP2;e;c g and fOIPe;c g will satisfy the continuous state inequality constraint (4). Theorem 2. Given 1; s 2 ½0; 1, for each e > 0, there exists a cðeÞ > 0 such that for all c; 0 < c < cðeÞ, any feasible solution of the model fOIP1;e;c ð1; sÞg; fOIP2;e;c g and fOIPe;c g is also a feasible solution of the model OIP1 ð1; sÞ; OIP2 and OIP. Remark 1. Theorem 2 ensures that the corresponding cðeÞ for each e in this sequence is finite. For constructing the optimization algorithm, when the parameter u1;s is not feasible, we can move the parameter towards e e;c ðl; u1;s ; 1; sÞ with respect to parameter u1;s . In this paper, the feasible region in the direction of the gradients of constraint G e e;c ðl; u1;s ; 1; sÞ. we develop a scheme for computing the gradients of constraint G e e;c ðl; u1;s ; 1; sÞ with Theorem 3. Given l 2 IN ; 1; s 2 ½0; 1, for each e > 0; c > 0, the derivatives of the constraint functionals G respect to the ith 2 I12 component of the parameter vector are

e e;c ðl; u1;s ; 1; sÞ @G :¼ @ui1;s

Z

tf

@Hðt; l; u1;s ; 1; s; v; xÞ dt; i ¼ 1; 2; . . . ; 12; @ui1;s

0

ð14Þ

where

Hðt; l; u1;s ; 1; s; v; xÞ :¼

16 Z X j¼1

0

tf

ue;j ðt; l; u1;s ; 1; sÞdt þ vT ðtÞhðxðt; xl0 ; u1;s ; 1; sÞ; u1;s ; 1; sÞ;

ð15Þ

and

vðtÞ :¼ ðv1 ðtÞ; v2 ðtÞ; . . . ; v8 ðtÞÞ;

ð16Þ

is the solution of the costate system (17)

v_ ðtÞ ¼ 

@Hðt; l; u1;s ; 1; s; v; xÞ ; @xðt; xl0 ; u1;s ; 1; sÞ

ð17Þ

with the boundary continuous

vð0Þ ¼ vðtf Þ ¼ ð0; 0; 0; 0; 0; 0; 0; 0ÞT :

ð18Þ

Proof. Let u1;s 2 U w ð1; sÞ be an arbitrary but fixed vector and di ; i 2 I12 , be an arbitrary real number. Define

ui;1;rs :¼ ðu11;s ; . . . ; ui1;s þ rdi ; . . . ; u12 1;s Þ; e e;c ðui;1;rs ; l; 1; sÞ can be where r > 0 is an arbitrarily small real number such that ui < ui1;s þ rdi < ui ; i 2 I12 . Therefore, G expressed as

2096

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

e e;c ðl; ui;r ; 1; sÞ :¼ c þ G 1;s

16 Z X j¼1

where

0

tf

ue;j ðt; l; ui;1;rs ; 1; sÞdt þ

Z 0

tf

_ xl0 ; ui;1;rs ; 1; sÞdt; vT ðtÞ½hðxðt; xl0 ; ui;1;rs ; 1; sÞ; ui;1;rs ; 1; sÞ  xðt;

ð19Þ

v is yet arbitrary. Thus, it follows that e e;c ðl; ui;1;rs ; 1; sÞ e e;c ðl; u1;s ; 1; sÞ dG @G jr¼0 ¼ di @ui1;s dr ( ) Z tf @Hðt; l; u1;s ; 1; s; v; xÞ @Hðt; l; u1;s ; 1; s; v; xÞ l l _ ¼ Dxðt; x0 ; u1;s ; 1; sÞ þ di  vðtÞDxðt; x0 ; u1;s ; 1; sÞ dt; @ui1;s @xðt; xl0 ; u1;s ; 1; sÞ 0

e e;c ðl; ui;r ; 1; sÞ :¼ DG 1;s

ð20Þ where Hðt; l; u1;s ; 1; s; v; xÞ is defined as in (15). Integrating (20) by parts and combining (15)–(19), we have

e e;c ðu1;s ; l; 1; sÞ @G di ¼ @ui1;s

Z

tf 0

@Hðt; l; u1;s ; 1; s; v; xÞ di dt: @ui1;s

ð21Þ

Since di is arbitrary, (14) follows readily from (21) and the proof is complete. h 5.3. Optimization algorithms Since the model OIP contains different parameter types, it would be hard to solve directly. To solve numerically OIP, various optimization routes, such as gradient-based algorithms [41], can be used. Nonetheless, all those techniques are only designed to find local optima. In order to overcome these difficulties, as mentioned in literature [32–34], particle swarm optimization (PSO) algorithm was originally developed for nonlinear optimization problems with different parameter types. Therefore, in this paper, PSO algorithm is applied to solve the model OIP. Firstly, we make a brief introduction to the basic PSO algorithm. Inspired by social behavior of bird flocking or fish schooling, in 1995, PSO algorithm was developed by Dr. Eberhart and Dr. Kennedy, which was a population on the basis of stochastic optimization technique [35]. Many similarities with evolutionary computation techniques were shared by PSO such as Genetic Algorithms (GA). A population of random solutions and search directions for optima by updating generations were used for initializing the system. However, unlike GA, PSO has no evolution operators such as crossover and mutation. In PSO, the potential solutions, called particles, fly through the problem space with directed velocity vectors to find better solutions. These velocity vectors have stochastic components and are dynamically adjusted based on historical and inter-particle information. Furthermore, previous studies [36] had indicated that the basic PSO algorithm converges quickly under all circumstances but slows its convergence speed down when approaching the optima, which exhibits notable prematurity. With this in mind, a novel particle swarm optimization algorithm based on particle migration (MPSO), which drew inspiration from the migratory behave in the nature when one particle travels from one swarm to another, was presented in detail in literature [36]. The proposed MPSO has the reputation of avoiding the local optima, and it can effectively take precautions against the premature convergence and significantly enhance diversity by information communication in the evolutionary process. The migratory scheme is given in Fig. 2. Since the parameters of different types are independent, we can split the OIP into a number of sub-problems. The main challenge is that each sub-problem contains a huge number of numerical computations of differential equations and

Fig. 2. The migration scheme [30].

2097

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

judgments of u1;s 2 U w ð1; sÞ. So, solving the whole OIP is intolerable for a serial program. In other words, it is impossible to work out within an acceptable time. In past several years, a parallel PSO algorithm had been widely applied in various fields [37–39]. In order to overcome above difficulties, the main steps of a parallel MPSO algorithm are as follows. 10    þ 1Þ þ j þ 1; i ¼ 0; 1; . . . ; m;  j ¼ 0; 1; . . . ; n  . Let Let hm ¼ 10  ¼ n  ; sk ¼ 0 þ jhn  ; k ¼ i  ðn  ; hn  ; m; n 2 Nþ ; 1k ¼ 0 þ ihm m  þ 1Þðn  þ 1Þ be the number of elements of D, and ð1k ; sk Þ 2 D be the kth (k=1; 2; . . . ; jDj) metabolic pathway. Denote jDj :¼ ðm CMNð1; sÞ as CMNðkÞ. Algorithm 1 Step 1: Allocate np processes for the algorithm. Let the population size be Psize . Then the number of particles is nw ¼ Psize =np on each process. Initialize data on the root process as follows.  Step 1.1: Read known data such as U ad ; W 0 . Set t f and a sufficient large number M sup .  Step 1.2: Read pathways set D, and let ð1k ; sk Þ 2 D as kth metabolic pathway, k ¼ 1; 2; . . . ; jDj.  Step 1.3: Read experimental data such as xl0 ; yl ðtÞ ¼ ðyl1 ðtÞ; yl2 ðtÞ; yl3 ðtÞ; yl4 ðtÞ; yl5 ðtÞÞ; t 2 N l (see Remark 2), l 2 IN . Step 2: Broadcast data on root process to all processes. Step 3: Initialize variables on the root process. opt  Step 3.1: Set J opt :¼ M sup ; k :¼ 0; k ¼ 1.  Step 3.2: Set Gopt :¼ M sup . Step 4: Execute the following procedure on each process, and denote the ID of process as p 2 Inp .  Step 4.1: Set Iter :¼ 1; Iv ar :¼ 1 and randomly generate nw particles with a uniform distribution on U ad . Denote the ter ter position and velocity of particles by uk;I and v k;I p;w p;w 2 ½V min ; V max ; w ¼ 1; 2; . . . ; nw , respectively. Here, 

u u  nv

 ; ½V min ; V max  :¼ ½ u nu v

  

is the range of velocities. Generally, nv :¼ 5 [39].

Step 4.2: Let M Iter be the maximum number of iterations, and M Iv ar be the maximum number where the optimal fitness does not improve continuously. Set P opt p;w :¼ M sup ; w 2 Inw . Step 4.3: Set w :¼ 1. ter Step 4.4: Put uk;I p;w in the OIP, solve the differential equations, defined by Eq. (2), by using the Euler algorithm and e e;c ðu1;s ; l; 1; sÞ by (13), 8l 2 IN , obtain the solution of CMNðkÞ, denoted by xð; xl ; uk;Iter ; kÞ; 8l 2 IN . Compute G 0



p;w

e e;c ðuk;Iter ; l; kÞ; SSDðu; 1; sÞ by (8), denoted by SSDðuk;Iter ; kÞ and U w ð1; sÞ by (9), denoted by U w ðkÞ for denoted by G p;w p;w kth metabolic pathway. e e;c ðuk;Iter ; l; kÞ P 0 and SSDðuk;Iter ; kÞ 2 U w ðkÞ, then compute Jðu; 1; sÞ by (10), denoted by Step 4.5: If 8l 2 IN ; G p;w p;w @e G ðu ;l;1;sÞ k;Iter k;Iter by Jðup;w ; kÞ. Otherwise, move the parameter up;w ðiÞ towards the feasible region in the direction of e;c@u1i ;s 1;s

@e G e;c ðup;wter ;l;kÞ k;I

(14), denoted by

k;I

@ðup;wter ðiÞÞ < P opt p;w ,

; 8i 2 I12 , by Theorem 3 with Armijo line searches, and go to Step 4.4.

opt opt k;Iter k;Iter l l k;Iter then set Popt  Step 4.6: If p;w :¼ Jðup;w ; kÞ; up;w :¼ up;w and xð; x0 ; up;w ; kÞ :¼ xð; x0 ; up;w ; kÞ; 8l 2 I N .  Step 4.7: If w < nw , then set w :¼ w þ 1, and go to Step 4.4. opt opt l Step 5: Gather P opt p;w ; up;w ; xð; x0 ; up;w ; kÞ; w ¼ 1; 2; . . . ; nw ; 8l 2 I N , from process p into the root process, denoted as c c c l c G ; u ; x ð; x0 ; u ; kÞ; c ¼ 1; 2; . . . ; P size ; 8l 2 IN . Step 6: Assign the root process to execute the tasks: If minc¼1;2;...;Psize Gc < Gopt then let c :¼ arg minc¼1;2;...;Psize Gc ; ter Jðuk;I p;w ; kÞ



Gopt :¼ Gc ; U opt :¼ U c



and set Iv ar :¼ 1. Otherwise, let c :¼ arg maxc¼1;2;...;Psize Gc ; :¼ G ; U :¼ U ; xwor;l :¼ xc ð; xl0 ; uc ; kÞ; 8l 2 IN . G Step 7: Broadcast U opt and Iv ar to all processes. Step 8: Migrating. If Iter %M I ¼¼ 0, then perform migration process Step 8.1 Step 8.3 on each process and denote the ID of the process as p 2 Inp ; Otherwise, go to Step 9.  Step 8.1: Assign a mutually exclusive index indicated as the target process number randomly for the pth 2 Inp process.  Step 8.2: Select mI migratory particles in the pth 2 Inp process, where mI ¼ nw r; r is the migration rate of the particle swarm, usually, r :¼ 0:5. The tournament selection is employed and the best particle should not be selected.  Step 8.3: If p ¼¼ np  1, then move the migratory particles from the pth 2 Inp process to the root process; otherwise, then move the migratory particles from the pth 2 Inp process to the pth þ 1 2 Inp process. Step 9: Execute the following procedure on each process.  Step 9.1: Update particles. Change the velocity and position of each particle according to the following rules: wor

c

wor



xopt;l :¼ xc ð; xl0 ; uc ; kÞ; 8l 2 IN

c

8 opt k;Iter ter > v k;Iter þ1 ðiÞ :¼ -ðIter Þv k;I > p;w ðiÞ þ c 1 r 1 ðup;w ðiÞ  up;w ðiÞÞ < p;w opt k;Iter i ¼ 1; . . . ; 12; w ¼ 1; . . . ; nw : þc2 r2 ðU ðiÞ  up;w ðiÞÞ; > > : uk;Iter þ1 ðiÞ :¼ uk;Iter ðiÞ þ v k;Iter þ1 ðiÞ; p;w

p;w

p;w

k;Iter ter where v k;I p;w ðiÞ and up;w ðiÞ represent the ith 2 I12 component of velocity and position of the wth 2 Inw particle of the pth 2 I np process at the Iter iteration, respectively. c1 ; c2 are the acceleration coefficients. r1 ; r 2 are the random numbers with the uniform distribution on [0,1]. In view of literature [36], it should be noted that the proper balance between global exploration

2098

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

and local exploitation is a key factor, which affects the performance of PSO algorithm, which is the proper balance between global exploration and local exploitation. By introducing the linearly decreasing inertia weight and time-varying acceleration coefficients into the original version of PSO, the performance of PSO has been remarkably improved. In solving OIP, after several times of numerical tests, we use the linearly decreasing inertia weight, which is given by Ma et al. [36]:

-ðIter Þ :¼ -max 

-max  -min M Iter

Iter ;

where -max and -min are the inertia weights of starting and ending, respectively. Usually the value of -ðIter Þ is varied between 0.9 and 0.4 [36]. Although the PSO algorithm with linearly decreasing inertia weight (LPSO) can locate satisfactory solution at a markedly fast speed, its ability to fine tune the optimum solution is limited. To cope with this problem, a novel strategy, which is implemented by changing the acceleration coefficients in such a manner that the cognitive component is reduced while the social component is increased as the search proceeds, is employed in literature [36]. The modification can be mathematically represented as follows:

c1 ¼ ðc1f  c1g Þ

Iter þ c1g ; MIter

c2 ¼ ðc2f  c2g Þ

Iter þ c2g ; MIter

where c1g ; c1f ; c2g and c2f are initial and final values of cognitive and social acceleration factors respectively, usually c1g ¼ c2f ¼ 2:5 and c1f ¼ c2g ¼ 0:5 [36]. When the parameters are out of their bounds, treated as follows:

( ter uk;I p;w ðiÞ

¼

ter ui ; if uk;I p;w ðiÞ < ui ;

 ter ui ; if uk;I p;w ðiÞ > ui ;

i ¼ 1; . . . ; 12:

 Step 9.2: If Iter < MIter and Iv ar < M Iv ar , then set Iter :¼ Iter þ 1; Iv ar :¼ Iv ar þ 1 and go to Step 4.3. Step 10: Execute the following procedure on the root process.  Step 10.1: Output Gopt and Gwor ; U opt and U wor ; xopt;l and xwor;l ; 8l 2 IN . If Gopt < J opt , then let J opt :¼ Gopt and opt k :¼ k.  Step 10.2: If k 6 jDj, then set k :¼ k þ 1 and go to Step 3.2. Otherwise, output J opt and 1kopt ; skopt , then stop. Table 1  n . The choice of parameters m;  m

10

20

30

40

 n mink2Ijðmþ1Þð Jðuk ; kÞ  nþ1Þj

10 0.006988

20 0.007012

30 0.007023

40 0.007029

1kopt skopt

0.9 0.9

0.85 0.90

0.90 0.87

0.875 0.9

0.4 0.35 0.3

k

J(u * ,k)

0.25 0.2 0.15 0.1 0.05 k=109

0 0

20

40

60

80

100

120

k Fig. 3. The values of the optimal index Jðuk ; kÞ; k ¼ 1; . . . ; jDj, where the marker is the optimal performance value of the model OIP the corresponding k ¼ 109.

2099

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

Remark 2. Let N l ¼ ftjt 2 ½0; t f andtis given moment in thelth 2 IN experimentg, where jN l j is the number of elements of the set N l . The formula (8) is integral forms, which would be difficult to compute directly. In addition, there are errors between yl ðtÞ; l 2 IN , and experiment data. Then, the formula (8) can be approximately calculated by

h i2 N X 3 X xj ðt; xl ; u; 1; sÞ  yl ðtÞ X 0 j 1 : SSDðu; 1; sÞ : ¼ PN : h i2 3 l¼1 jNl j l¼1 j¼1 t2Nl yl ðtÞ

ð22Þ

j

The formula (10) is integral forms, which would be difficult to compute directly. To overcome these difficulties, our methods are transformation of numerical integration and to divide the perturbation space U with uniformity into g subdivisions. Then, by applying Mean-Value Theorem on every subdivision, we can find duðiÞ such that (10) can be approximately calculated by

Jðu; 1; sÞ : ¼_

tf  g 

1 XN l¼1

g N X X

jNl j

P

t2N l

kxðt;xl0 ;uþduðiÞ ;1;sÞxðt;xl0 ;u;1;sÞk kxðt;xl0 ;u;1;sÞk

ð23Þ

:

kduðiÞ k

l¼1 i¼1

Table 2 ~ ; u109 . The value of parameters u u1 109

u2 109

u3 109

u4 109

u5 109

u6 109

34.0297 u7 109 1.98077 ~1 u 30 ~7 u 1.5

107.204 u8 109 17.7808 ~2 u 100 ~8 u 18

150.417 u9 109 10.6496 ~3 u 170 ~9 u 12

4.71128 u10 109 5263.6 ~4 u 4 ~ 10 u 5000

115.59 u11 109 44.8808 ~5 u 120 ~ 11 u 40

57.2065 u12 109 68.8345 ~6 u 60 ~ 12 u 70

4

2

1

2

3

4

5

6

600

400

200

0

7

intracellular glycerol(mmol/L)

intracellular glycerol(mmol/L)

biomass(g/L)

3

0

800

800 computational results experimental data

0

2

time(h)

400 300 200 100 0

3

4

5

6

300

250

250

200 150 100

0

7

0

2

4 time(h)

6

200 150 100 50 2

3

4

5 time(h)

6

7

4 time(h)

6

8

0

2

4 time(h)

6

8

0

2

4 time(h)

6

8

150 100

250

200 150 100 50 0

2

200

0

8

intracellular 1,3−PD(mmol/L)

250

intracellular 1,3−PD(mmol/L)

extracelluar 1,3−PD(mmol/L)

computational results experimental data

0

50

250

300

200

0

8

50 2

400

300

time(h)

0

6

3−HPA(mmol/L)

computational results experimental data

3−HPA(mmol/L)

extracelluar glycerol(mmol/L)

500

4 time(h)

600

0

2

4 time(h)

6

8

200 150 100 50 0

Fig. 4. The results of the 109th metabolic pathway of initial glycerol concentration of 441.337 mmol=L and biomass concentration of 0.2025 g=L.

2100

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

The formula (23) represents the average rate of change with respect to some perturbations duðiÞ . Considering the amount of calculation, we take g ¼ 10. In other words, some perturbations are carried out with each iteration, and the average rate of change is taken under the perturbations as the approximation of Jðu; 1; sÞ. Combining Algorithm 1 with Theorem 2, we can develop the following algorithm to solve the (OIP). Algorithm 2 Step 1: Choose initial values of

e > 0; c > 0.

Step 2: Solve Problem ðOIP e;c Þ using Algorithm 1 to give ukopt ;e;c and Step 3: Check feasibility of

g j ðt; l; ukopt ;e;c ;

1kopt ; skopt .

1kopt ; skopt Þ for 8t 2 ½0; tf ; j 2 I16 .

, where c  is a prespecified positive constant, go to Step Step 4: If ukopt ;e;c is feasible, go to Step 5. Otherwise, set c :¼ ac. If c > c 2. Otherwise go to Step 6. Step 5: Set e :¼ be. If e > e, where e is a prespecified positive constant, go to Step 2. Otherwise go to Step 6. Step 6: Output ukopt ;e;c ; 1kopt ; skopt and stop. Remark 3. In Algorithm 2, e is a parameter controlling the accuracy of the smoothing approximation. c is a parameter con are two predefined parameters ensuring the termination of the algorithm. trolling the feasibility of the constraint(4). e and c . Especially, the parameters a Remark 4. It is important for the validity of Algorithm 2 to choose the parameters a; b; e and c  are two sufficiently small values such that Algorithm 2 can be terminated. and b must be chosen less than 1. e and c 6. Numerical results The below parameters are derived empirically after numerous experiments: tf :¼ 7:75; n :¼ 0:3; N :¼ 3; N 1 :¼ f2; 3; 4; 5; 6; 6:5g; N 2 :¼ f1; 2; 3; 4; 5; 6; 6:67g; N 3 :¼ f1:5; 2:5; 3:5; 4:75; 5:75; 6:75; 7:75g; x10 :¼ f0:2025; 441:337; 0; 0; 0; 0; 0; 0g; x20 :¼ f0:173; 402:9348; 0; 0; 0; 0; 0; 0g; x30 :¼ f0:2245; 509:8913; 0; 0; 0; 0; 0; 0g; The system CMNðkÞ; k ¼ 1; 2; . . . ; jDj, is solved

2

1

0

1

2

3

4 time(h)

5

6

3−HPA(mmol/L)

extracelluar glycerol(mmol/L)

400

200

100

0

computational results experimental data 1

2

3

4 time(h)

0

2

4 time(h)

6

5

6

7

computational results experimental data

200 150 100 50 1

2

3

4 time(h)

5

6

7

400

200

0

8

300

250

250

200 150 100

0

2

4 time(h)

6

2

4 time(h)

6

8

0

2

4 time(h)

6

8

0

2

4 time(h)

6

8

200 150 100

0

8

250

200 150 100 50 0

0

50

250

250

600

300

0

intracellular 1,3−PD(mmol/L)

extracelluar 1,3−PD(mmol/L)

200

50

300

0

400

0

7

300

600

intracellular 1,3−PD(mmol/L)

biomass(g/L)

3

800

3−HPA(mmol/L)

computational results experimental data

intracellular glycerol(mmol/L)

800

intracellular glycerol(mmol/L)

4

0

2

4 time(h)

6

8

200 150 100 50 0

Fig. 5. The results of the 109th metabolic pathway of initial glycerol concentration of 402.9348 mmol=L and biomass concentration of 0.173 g=L.

2101

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103 5

3 2 1

2

4 time(h)

6

600 400 200

0

2

4 time(h)

6

800 600 400 200 0

8

300

300

400

250

250

300 200 100

computational results experimental data 0

2

4 time(h)

200 150 100 50

6

0

8

0

2

4 time(h)

6

400 300 200 100

0

2

4 time(h)

6

8

4 time(h)

6

8

0

2

4 time(h)

6

8

0

2

4 time(h)

6

8

150 100

400

300

200

100

0

2

200

0

8

intracellular 1,3−PD(mmol/L)

intracellular 1,3−PD(mmol/L)

computational results experimental data

0

50

400

500

0

3−HPA(mmol/L)

500

0

extracelluar 1,3−PD(mmol/L)

800

0

8

3−HPA(mmol/L)

extracelluar glycerol(mmol/L)

0

intracellular glycerol(mmol/L)

intracellular glycerol(mmol/L)

biomass(g/L)

4

0

1000

1000 computational results experimental data

0

2

4 time(h)

6

8

300

200

100

0

Fig. 6. The results of the 109th metabolic pathway of initial glycerol concentration of 509.8913 mmol=L and biomass concentration of 0.2245 g=L.

by using Euler method with a step size of 1/3600 (h). In Algorithm 1, M I :¼ 5; M Iv ar :¼ 200; M Iter :¼ 600; Psize :¼ 80; np :¼ 40.  :¼ 1:0 107 . In Algorithm 2, e :¼ 0:1; c :¼ 0:01; a :¼ 0:1; b :¼ 0:01; e :¼ 1:0 107 ; c 12  Let uk ¼ ðu1k ; u2k ; . . . ; u12 Þ 2 R be system parameters and u be the optimal parameter of the kth ; k ¼ 1; 2; . . . ; jDj metak k bolic pathway. Denote Jðu; 1; sÞ by (10) as Jðuk ; kÞ and U w ð1; sÞ by (9) as U w ðkÞ for kth, k ¼ 1; 2; . . . ; jDj, metabolic pathway.  ¼n  ¼ 10. From above After numerous experiments, it follows from Table 1 and the amount of calculation that we take m  þ 1Þðn  þ 1Þ ¼ 121 metabolic pathways, 1452 ¼ jDj 12 system parameters and description, we can see that jDj ¼ ðm 484 ¼ jDj 4 pathway parameters are involved in the model OIP. The model OIP is solved by using Algorithms 1 and 2, including 9979188 numerical computations of differential equations on Lenovo DeepComp 1800 PC-cluster Server composed of 16 nodes. Each node is equipped with two Intel 5420 CPU (4 Core, 64-bit, clocked at 2.5 GHz) and 8 GB memory. By employing 40 processors, the total computational time was about 17.75 h. The best biological robustness of the kth ; k ¼ 1; 2; . . . ; jDj metabolic pathway are shown in Fig. 3 by Algorithms 1 and 2. Form Fig. 3, for the 1st, 56th, 111th, 116th and 121st metabolic pathways, their biological robustness are set to be þ1 because the sets U w ð1Þ; U w ð56Þ; U w ð111Þ; U w ð116Þ; U w ð121Þ, are respectively empty even if the n is large enough, which indicates that the 1st, 56th, 111th, 116th and 121st metabolic pathways, are all very poor in representing the actual fermentation process. That is, it follows from U w ð1Þ; U w ð56Þ; U w ð111Þ that the transport way of 1,3-PD cannot be active transport. It follows from U w ð111Þ; U w ð116Þ; U w ð121Þ that the transport way of glycerol cannot be passive diffusion. So opt we conclude that the metabolic pathway k ¼ 109 is the most possible. That is, it is most possible that glycerol and 1,3PD passes the cell membrane by active transport coupled with passive diffusion (i.e. 1kopt ¼ 0:9; 1  1kopt ¼ 0:1; skopt ¼ 0:9; 1  skopt ¼ 0:1). The identified system parameters of the k ¼ 109 metabolic pathway are shown in Table 2 and in Figs. 4–6. ~ ; u109 ðwhere any u ~ – u109 Þ, are given in Table 2. As you can be seen in Figs. 4–6, on the left of three pictures, the dotted line u denotes computational results, and the solidline curve denotes experimental data. In Figs. 4–6, the middle and the right of three pictures denote simulated results of intracellular glycerol, intracellular 3-HPA and intracellular 1,3-PD concentrations ~ ; 0:02Þ ðwhere any u ~ – u109 Þ and of the dynamical system (1) with g ¼ 1000 the uniformly distributed parameters in Bðu  Bðu109 ; 0:02Þ respectively, at different initial glycerol and biomass concentration of 441.337 mmol/L and 0.2025 g/L, 402.9348 mmol/L and 0.173 g/L, together with 509.8913 mmol/L and 0.2245 g/L in order. The white denotes concentra-

2102

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

tions of intracellular substances when uk – u109 and uk ¼ u109 , respectively; and the black point line denotes concentrations of intracellular substances when u109 is disturbed in the vicinity of uk – u109 and uk ¼ u109 , respectively. We can see from Figs. 4–6 that system (1) is less sensitive to fluctuations on parameters uk ¼ u109 . Based on Definition 2, when k ¼ 109, ~ ; 1; sÞ; 81; s 2 f0; 0:1; 0:2; 0:3; 0:4; 0:5; 0:6; 0:7; 0:8; 0:9; 1g. Numerical intei.e. 1 ¼ 0:9; s ¼ 0:9, we have Jðu109 ; 0:9; 0:9Þ 6 Jðu ~ 2 U ad , which grations illustrate that the problem with parameter vector u109 is more robust than any parameter vector u is revealed from Figs. 4–6. 7. Discussions and conclusions In this paper, we establish a complex metabolic system, in which 121 possible metabolic pathways of glycerol and 1,3-PD are involved, to describe the batch fermentation of glycerol by K. pneumoniae. We then discuss some basic properties of the metabolic system. Inspired by the qualitative description of the biological robustness given by Kitano and h-stability, we give a quantitative definition of biological robustness of intracellular substances, which is different from quantitative description of biological robustness of continue culture in the approximately stable state, for the overall process of batch culture. Taking biological robustness of intracellular substances as a cost function, we present an identification model, in which 1452 system parameters and 484 pathway parameters are involved, subject to the metabolic system, continuous state constraints and parameter constraints. Subsequently, its identifiability is proved. The constraint transcription and smoothing techniques are applied to deal with the continuous state constraints in the identification model. To determine the most possible metabolic pathway, we construct a parallel optimization algorithm, based on particle swarm optimization migration (MPSO), for solving the identification model. Through numerical calculation, we conclude that it is most possible that both glycerol and 1,3-PD pass the cell membrane by active transport and passive diffusion with different weights. In a future work, our efforts will focus on multistage analysis and optimal control of glycerol input in batch culture. Additionally, necessary condition of optimality for the system identification model is waiting for publication in another paper and the stability of Algorithm 1 will be discussed at full length in the future. Acknowledgments This work was supported by 863 Program (Grant No. 2007AA02Z208), 973 Program (Grant No. 2007CB714304), the Fundamental Research Funds for the Central Universities (Grant No. DUTTX2011106), the National Natural Science Foundation of China (Grant No. 11171050), the National Natural Science Foundation for the Youth of China (Grant No. 11301051). References [1] Witt U, Miiller RJ, Augusta J, Widdecke H, Deckwer WD. Synthesis, properties and biodegradability of polyesters based on 1,3-propanediol. Macromol Chem Phys 1994;195:793–802. [2] Biebl H, Marten S, Hippe H, Deckwer WD. Glycerol conversion to 1,3-propanediol by newly isolated clostridia. Appl Microbiol Biotechnol 1992;36:592–7. [3] Nakamura CE, Whited GM. Metabolic engineering for the microbial production of 1,3-propanediol. Curr Opin Biotechnol 2003;14:454–9. [4] Gtinzel B. Mikrobielle herstellung von 1,3-propandiol durch Clostridium butyricum und adsorptive Abtremutng von Diolen. Ph.D. dissertation, TU Braunschweig, Germany; 1991. [5] Biebl B, Menzel K, Zeng AP, Deckwer WD. Microbial production of 1,3-propanediol. Appl Microbiol Biotech 1999;52:297–8. [6] Zeng AP, Rose A, Biebl H, Tag C, Guenzel B, Deckwer WD. Multiple product inhibition and growth modeling of Clostridium butyricum and Klebsiella pneumoniae in fermentation. Biotechnol Bioeng 1994;44:902–11. [7] Zeng AP. A kinetic model for product formation of microbial and mammalian cells. Biotechnol Bioeng 1995;46:314–24. [8] Wang L, Xiu ZL, Zhang YD, Feng EM. Optimal control for multistage nonlinear dynamic system of microbial bioconversion in batch culture. J Appl Math 2011. http://dx.doi.org/10.1155/2011/624516. 11p [Article ID 624516]. [9] Wang L, Feng EM, Ye JX, Xiu ZL. An improved model for multistage simulation of glycerol fermentation in batch culture and its parameter identification. Nonlinear Anal Hybrid Syst 2009;3:455–62. [10] Wang J, Ye JX, Yin HC, Feng EM, Wang L. Sensitivity analysis and identification of kinetic parameters in batch fermentation of glycerol. J Comput Appl Math 2012;236:2268–76. [11] Wang L, Feng EM, Xiu ZL. Modeling nonlinear stochastic kinetic system and stochastic optimal control of microbial bioconversion process in batch culture. Nonlinear Anal Model 2013;18:99–111. [12] Barkai N, Leibler S. Robustness in simple biochemical networks. Nature 1997;387:913–7. [13] Kitano H. Biological robustness. Nat Rev Genet 2004;5:826–37. [14] Kitano H. Towards a theory of biological robustness. Mol Syst Biol 2007;3:137. [15] Marhl M, Perc M. Determining the flexibility of regular and chaotic attractors. Chaos Soliton Fractals 2006;28:822–33. [16] Perc M, Marhl M. Sensitivity and flexibility of regular and chaotic calcium oscillations. Biophys Chem 2003;104:509–22. [17] Ozer M, Uzuntarla M, Perc M, Graham LJ. Spike latency and jitter of neuronal membrane patches with stochastic Hodgkin–Huxley channels. J Theor Biol 2009;261:83–92. [18] Perc M, Marhl M. Noise enhances robustness of intracellular Ca2+ oscillations. Phys Lett A 2003;316:304–10. [19] Marhl M, Gosak M, Perc M, Roux E. Importance of cell variability for calcium signaling in rat airway myocytes. Biophys Chem 2010;148:42–50. [20] Stoz˘er A, Gosak M, Dolens˘ek J, Perc M, Marhl M, Rupnik MS, Koros˘ak D. Functional connectivity in islets of Langerhans from mouse pancreas tissue slices. Comput Biol 2013;9:e1002923. [21] Bodenstein C, Gosak M, Schuster S, Marhl M, Perc M. Modeling the seasonal adaptation of circadian clocks by changes in the network structure of the suprachiasmatic nucleus. Comput Biol 2012;8:e1002697 [PLoS ONE]. [22] Zhang YD, Feng EM, Xiu ZL. Robust analysis of hybrid dynamical systems for 1,3-propanediol transport mechanisms in microbial continuous fermentation. Math Comput Model 2011;54:3164–71. [23] Yan HH, Zhang X, Ye JX, Feng EM. Identification and robustness analysis of nonlinear hybrid dynamical system concerning glycerol transport mechanism. Comput Chem Eng 2012;40:171–80.

J. Yuan et al. / Commun Nonlinear Sci Numer Simulat 19 (2014) 2088–2103

2103

[24] Sun YQ, Qi WT, Teng H, Xiu ZL, Zeng AP. Mathematica modeling of glycerol fermentation by Klebsiella pneumoniae: concerning enzyme-catalytic reductive pathway and transport of glycerol and 1,3-propanediol across cell membrane. Biochem Eng J 2008;38:22–32. [25] Sun J, Heuvel J, Soucaille P, Qu Y, Zeng AP. Comparative genomic analysis of DHA regulon and related genes for anaerobic glycerol metabolism in bacteria. Biotechnol Prog 2003;19:263–72. [26] Choi SK, Koo NJ. H-stability for nonlinear perturbed systems. Ann Differ Equ 1995;11:1–9. [27] Teo KL, Goh CJ, Wong KH. A unified computational approach to optimal control problems. Essex: Long Scientific Technical; 1991. [28] Torres NV, Voit EO. Pathway analysis and optimization in metabolic engineering. Cambridge: Cambridge University Press; 2002. [29] Jiang ZG, Yuan JL, Feng EM. Robust identification and its properties of nonlinear bilevel multi-stage dynamic system. Appl Math Comput 2013;219:6979–85. [30] Wang L. Determining the transport mechanism of an enzyme-catalytic complex metabolic network based on biological robustness. Bioprocess Biosyst Eng 2013;36:433–41. [31] Wang L. Modelling and regularity of nonlinear impulsive switching dynamical system in fed-batch culture. Abstr Appl Anal 2012;2012. http:// dx.doi.org/10.1155/2012/295627. 15p [Article ID 295627]. [32] Gao WF, Liu SY, Huang LL. Particle swarm optimization with chaotic opposition-based population initialization and stochastic search technique. Commun Nonlinear Sci Numer Simul 2012;17:4316–27. [33] Gandomi AH, Yun GJ, Yang XS, Talatahari S. Chaos-enhanced accelerated particle swarm optimization. Commun Nonlinear Sci Numer Simul 2013;18:327–40. [34] Zhang J, Zhang C, Chu T, Perc M. Resolution of the stochastic strategy spatial prisoner’s dilemma by means of particle swarm optimization. Comput Biol 2011;6(7):e21787 [PLoS ONE]. [35] Kennedy J, Eberhart R. Particle swarm optimization. In: Proceedings of IEEE international conference on neural networks. Piscataway, NJ; 1995. p. 1942–8. [36] Ma G, Zhou W, Chang XL. A novel particle swarm optimization algorithm based on particle migration. Appl Math Comput 2012;218:6620–6. [37] Koh B, George AD, Haftka RT, Fregly BJ. Parallel asynchronous particle swarm optimization. Int J Numer Methods Eng 2006;76:578–95. [38] Wang Y, Wong KW, Xiao D. Parallel hash function construction based on coupled map lattices. Commun Nonlinear Sci Numer Simul 2011;16:2810–21. [39] Zhang YL, Gallipoli D, Augarde C. Parameter identification for elasto-plastic modelling of unsaturated soils from pressuremeter tests by parallel modified particle swarm optimization. Comput Geotech 2013;48:293–303. [40] Goberna MA, López MA. Semi-infinite programming recent advances. Dordrecht: Kluwer; 2001. [41] Polak E. Optimization algorithms and consistent approximations. New York: Springer-Verlag; 1997.