Identification and robustness analysis of nonlinear hybrid dynamical system of genetic regulation in continuous culture

In this paper, we present a framework to infer the possible transmembrane transport of intracellular substances. Considering four key enzymes, a modified fourteen-dimensional nonlinear hybrid dynamic system is established to describe the microbial continuous culture with enzyme-catalytic and genetic regulation. A novel quantitative definition of biological robustness is proposed to characterize the system's resilience when system parameters were perturbed. It not only considers the expectation of system output data after parameter disturbance but also considers the influence of the variance of these data. In this way, the definition can be used as an objective function of the system identification model due to the lack of data on the concentration of intracellular substances. Then, we design a parallel computing method to solve the system identification model. Numerical results indicate that the most likely transmembrane mode of transport is active transport coupling with passive diffusion for glycerol and 1, 3-propanediol.


Zhilong Xiu
School of Life Science and Biotechnology, Dalian University of Technology Dalian, Liaoning 116024, China

(Communicated by Changzhi Wu)
Abstract. In this paper, we present a framework to infer the possible transmembrane transport of intracellular substances. Considering four key enzymes, a modified fourteen-dimensional nonlinear hybrid dynamic system is established to describe the microbial continuous culture with enzyme-catalytic and genetic regulation. A novel quantitative definition of biological robustness is proposed to characterize the system's resilience when system parameters were perturbed. It not only considers the expectation of system output data after parameter disturbance but also considers the influence of the variance of these data. In this way, the definition can be used as an objective function of the system identification model due to the lack of data on the concentration of intracellular substances. Then, we design a parallel computing method to solve the system identification model. Numerical results indicate that the most likely transmembrane mode of transport is active transport coupling with passive diffusion for glycerol and 1,3-propanediol.
1. Introduction. The bioconversion of glycerol to 1,3-propanediol (1,3-PD) has a wide range of potential applications on a large commercial scale, especially as a monomer for polyesters, polyethers and polyurethanes, due to its lower cost, higher production rate and no pollution [2,35]. There are three typical cultures for microbial fermentation of glycerol, including batch, continuous and fed-batch cultures. Over the past years, a number of models have been proposed to describe the bio-dissimilation process of glycerol to 1,3-PD by Klebsiella pneumoniae (K. pneumoniae).
In recent years, some nonlinear dynamical systems have been proposed to describe the microbial fermentation process mentioned above. A well-cited suitable dynamic model was first proposed by Zeng et al. [43]. In this model, only the concentrations of microorganisms, glycerol, 1,3-PD, and two other major by-products were considered. Then, Xiu et al. [36,37] made further improvements based on this model to analyze some qualitative and quantitative characteristics of the modified model. Based on this improved mathematical model, some scholars have studied the parameter identification problem, stability of equilibria, phenomena of oscillation and hysteresis for the continuous fermentation [5,6,14,16,17,18,21,22,39]. However, the above literatures didn't consider the effect of changes in the concentration of intracellular substances on the microbial fermentation reaction. With the deepening of research, such as the understanding of intracellular transport of substances across the membrane can be very helpful for the research in genetic engineering and other related aspects, the impact of intracellular substances cannot be ignored. Therefore, in 2008, Sun et al. [27] proposed a novel mathematical model, in which the concentration changes of both extracellular substances and intracellular substances were taken into consideration. Thereafter, a large number of further investigations have been carried out, including sensitivity analysis [31], pathway identification problems [30,38,44], stochastic optimal control [33] and robust parameter identification [42].
In 2012, Sun et al. [28] proposed a fourteen-dimensional nonlinear dynamic system concerning the repression of the dha regulon and two key enzymes by 3hydroxypropionaldehy (3-HPA) to describe the continuous culture, which is more accurate and complex than eight-dimensional model. On the basis of the fourteendimensional model, Yuan et al. [41] establish a nonlinear hybrid dynamical system with uncertain system parameters and pathway parameters for describing the process of batch culture. Gao et al. [7] propose a novel global sensitivity analysis method to reduce the parameters in the parameter identification model. Considering the three possible transport mechanisms of 1,3-PD across cell membrane. The transmembrane transport of intracellular glycerol was not considered in the above literature, Guo et al. [8] establish a nonlinear hybrid dynamical system about continuous culture, as this would impose a huge computational burden on the subsequent analysis.
Since it is difficult to accurately measure the concentration of intracellular substances, it is inevitable to determine the most reasonable glycerol metabolism mechanism by means of numerical analysis. Hiroaki [12,13] proposed a popular concept of biological robustness to judge the reliability of numerical solution for the concentrations of intracellular substances. Biological robustness is a basic characterization of biological systems, which has also been accepted by researchers in the field [19,24,25,34]. Many researchers have tried to give a quantitative definition of biological robustness and achieved some desired results [32,40]. The core idea of quantitative definition is to quantify the changes in system performance by perturbing system parameters. In this paper, we establish a parameter identification model with a novel quantitative definition of biological robustness. And through this identification model to infer the most likely glycerol transmembrane transport. The novel quantitative definition not only considers the expectation of system output data after parameter disturbance but also considers the influence of the variance of these data. In this way, the definition can be used as an objective function of the system identification model due to the lack of data on the concentration of intracellular substances, the relative error of extracellular substances and the approximate stability of dynamical systems as constraint conditions. Finally, we design a parallel computing method to solve the system identification model. Then, we obtain the optimal pathway and parameters of the corresponding system. This paper is organized as follows. In Section 2, a nonlinear dynamical system of enzyme-catalytic and genetic regulation of continuous culture is presented. In Section 3, some important properties of the system are discussed. In Section 4, a new quantitative description of biological robustness is presented and with this presented biological robustness as performance index, a system identification model is established. Section 5 constructs a parallel particles warm optimization algorithm for the system identification model. In Section 6, a numerical result is given to infer the parameter identification model. Conclusions are presented at the end of the paper.
2. Nonlinear hybrid dynamical system with genetic regulation. The biochemical process under consideration is the dissociation of glycerol in Klebsiella pneumoniae, which is a disproportionation process that occurs through coupled oxidation and reduction pathways, as shown in Fig.1. During glycerol metabolism, glycerol is transformed from the extracellular environment to the intracellular environment and then converted to 3-hydroxypropionaldehyde (3-HPA) by glycerol dehydratase (GDHt). The important intermediate 3-HPA was further reduced to 1,3-PD by enzyme 1,3-propanediol oxidoreductase (PDOR). Finally, the by-products of the target product 1,3-PD and coupled oxidation are transformed from the intracellular environment to the extracellular environment. Based on the above process, we propose a fourteen-dimensional nonlinear dynamic system with two key reductases to describe continuous culture [28].
Since the transport mechanisms of glycerol and 1,3-PD across cell membrane in the reductive pathway are still unclear, we shall regard the metabolic pathway with a combination of possible transport mechanisms of glycerol and 1,3-PD as a possible metabolic system, which is convenient to be represented by a simple metabolic network. Taking all possible transport mechanisms of glycerol and 1,3-PD into consideration, we can obtain nine possible combinations of transports mechanisms, which correspond to nine possible simple metabolic networks (or nine possible metabolic systems), denoted by NHDS(l 1 ), NHDS(l 2 ), · · · , NHDS(l 9 ), see Table 1.
According to the experiment process, we make the following assumptions: (H1) The concentrations of reactants are uniform in reactor, time delay and nonuniform space distribution are ignored.
respectively denote the concentration of biomass, extracellular glycerol, extracellular 1,3-PD, acetic, ethanol, intracellular glycerol, intracellular 3-HPA, intracellular 1,3-PD, mR, R, mGDHt, GDHt, mPDOR and PDOR at time t, respectively. To simplify notations, I n := {1, 2, · · · , n} and x i := x i (t). Let N be the total number of experiments carried out under different dilution rates and initial glycerol concentrations (D j , C s0,j ), j ∈ I N , and p j,k := p 1 j,k , p 2 j,k , . . ., p 31 j,k T , k ∈ I 9 be the kinetic parameter vector to be determined. Under the assumptions (H1) and (H2), the continuous fermentation process of glycerol by K. pneumoniae of the kth metabolic system in the jth experiment, denoted by NHDS(j, l k ), can be described as follows: where D j and C s0,j are respectively the dilution rate and the glycerol concentration in the feed medium in the jth experiment. The specific cellular growth rate µ in Eq.(1) can be given as [43]: For K. pneumoniae grown under anaerobic conditions at 37 • C and pH value 7.0, the maximum specific growth rate µ m and Monod constant K s in Eq.(13) take the values of 0.67h −1 and 0.28mmol/L , where, K mG and K mP denote Michaelis-Menten constants of enzymes GDHt and enzymes PDOR, which are 0.53mmol/L and 0.14mmol/L respectively. The transport system of glycerol and 1,3-PD are considered, according to Monod equation and Fick diffusion law, the specific consumption rate of glycerol and specific product formation of 1,3-PD can be expressed by where l i k ∈ {0, 1}, i ∈ I 4 , k ∈ I 9 , l 1 k = 0 (l 2 k = 0) and l 3 k = 0 (l 4 k = 0) indicates glycerol and 1,3-PD pass the membrane the passive diffusion (resp. active transport) doesn't exist, whereas l 1 k = 1 (l 2 k = 1) and l 3 k = 1 (l 4 k = 1) represents the

QI YANG, LEI WANG, ENMIN FENG, HONGCHAO YIN AND ZHILONG XIU
presence of active transport of glycerol and 1,3-PD (resp. passive diffusion), respectively. According to the actual fermentation process, we have max{l 1 k , l 2 k } = 1 and max{l 3 k , l 4 k } = 1. Let l k := (l 1 k , l 2 k , l 3 k , l 4 k ) T ∈ S l := {l 1 , l 2 , · · · , l 9 }, k ∈ I 9 , then l k corresponds to the kth metabolic system NHDS(l k ) in the complex metabolic network NHDS. For details, see Table 1. Table 1. Transport mechanisms of glycerol and 1, 3-PD of metabolic system N HDS(l k ), corresponding to parameter vector l k , k ∈ I 9 . Abbreviations: A, active transport; P, passive diffusion; AP, passive diffusion coupled with active transport.
In addition, the indicator function defined as The uptake of extracellular glycerol is considered as "black box" model, its specific consumption rate expressed by [27]: The specific formation rate of extracellular acetate and ethanol can be expressed respectively as follows [27]: The parameters in Eqs.(16)- (18) were estimated by the previous work [37], as shown in Table 2. To simplify notation, we denote p j,k as p, the nonlinear hybrid dynamical system NHDS(j, l k ), can be described by the combination of the following two systems: where Denote the time intervals of batch and continuous periods by The system (19) describes the batch fermentation process, in which x 0 is the initial state, x b is the terminal state, and the dilution rates to be zero. The system (20) describes the continuous fermentation process starting from the terminal state of the system (19) with prescribed dilution rate D j = 0 . According to the practical culture process, the concentrations of biomass, substrate, intermediates and products will not exceed their critical concentrations. Let x * i and x * i be the lower and upper bounds of the state vector x i , so we consider the properties of the system on: x * i = (0.001, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) T , Let p * and p * be the lower and upper bounds of the parameter vector p j,k , denote P ad be the allowable set of vector p j,k , which can be given by [28]: where | p i,0 | := 0.5p i,0 , 3. Properties of the nonlinear hybrid dynamical system. The existence, uniqueness and continuous differentiability of the solution to the system NHDS(j, l k ) with respect to parameters are explored in this section. In addition, we denote Euclidean norm by · in this paper.
According to the factual fermentation mechanism, we have the following assumptions: (H3) The set W ad , D c and P ad are all nonempty, bounded and closed convex sets.
(H4) The absolute differences of concentrations between extracellular and intracellular glycerol concentration and that of 1, 3-PD concentration are uniformly bounded, i.e., ∃M 1 > 0, M 2 > 0, such that Proof. The conclusions of the property can be easily obtained by using the methods similar to the proof of Property 1 in [7].
Under the assumptions (H3) and (H4) and Property 1, we can obtain the following properties.
Property 2. For given (j, l k ) ∈ I N × S l and p∈P ad , there exists a unique solution to the system NHDS(j, l k ), denoted by x(·; x 0 , p, j, l k ). Furthermore, x(·; x 0 , p, j, l k ) is continuous in p on P ad .
For given (j, l k ) ∈ I N × S l , p∈P ad , x 0 ∈ W ad , we define the following three sets: 4. Biological robustness analysis and parameter identification model. The purpose of this section is to establish a reasonable criterion to infer the most impossible metabolic mechanism. In continuous fermentation experiments, although the steady-state data of extracellular substances can be measured, the intracellular metabolic mechanisms are not completely clear, and it is difficult to measure the concentration of the intracellular substances and the accuracy of measurement.
To deal with this problem, a novel quantitative definition of biological robustness is proposed to characterize the system's resilience when system parameters were perturbed. It not only considers the expectation of system output data after parameter disturbance but also considers the influence of the variance of these data. Thus, this quantitative definition can be used as part of the objective function in the system identification model to characterize the lack of data for intracellular concentrations. Then we discuss the existence of the optimal solution to the identification problem. Some basic definitions are given as follows.  [31]) For given (j, l k ) ∈ I N × S l , a state vector x(t s ; x 0 , p, j, l k ) ∈ W ad is said to be an approximately stable solution within a precision η > 0, if there exists p ∈ D p (j, l k ) and a t s ∈ [t 0 , t f ] such that x(t; x 0 , p, j, l k ) is a solution of system NHDS(j, l k ) satisfying where t s := inf t σ : f i (t; x, p, j, l k ) ≤ η, ∀ t∈ [t σ , t f ] , t s is the steady state time.
Definition 4.2. (Relative Error of Extracellular Substances) Given (j, l k ) ∈ I N ×S l and p ∈ D p (j, l k ), let x i (t s ; x 0 , p, j, l k ), i ∈ I 3 be the ith component of approximately stable state by the system NHDS(j, l k ) and y i (j)∈R 3 , i ∈ I 3 be the j th experimental extracellular concentration (i.e., biomass , glycerol and 1,3-PD) measured at the steady stage. Then, the relative error of extracellular substances is defined as: .
The set of feasible parameter vector is defined as follows: where the ε 0 is a given small constant, D ps (j, l k ) ⊆ D p (j, l k ) is nonempty and compact in R 31 . For givenp := (p 1 j,k ,p 2 j,k , · · · ,p 31 j,k ) T ∈ D ps (j, l k ) and p := (p 1 j,k , p 2 j,k , · · · , p 31 j,k ) T ∈ D ps (j, l k ). We define the rectanglar field of dimension normalization, denoted by B(p, δ), can be expressed as follows: where δ > 0 is sufficiently small to measure the magnitude of the parameter disturbances. Obviously, B(p, δ) is compact in R 31 .
Definition 4.3. (Relative deviation with respect to parameters perturbation) Given (j, l k ) ∈ I N ×S l , assume that the system reaches an approximately stable state at t 1 s and t 2 s with the kinetic parameter vector taking values p , p ∈ B(p, δ), respectively. Let t s := max{t 1 s , t 2 s } denote the corresponding approximately stable solutions by x i (t s ; p , j, l k ) and x i (t s ; p , j, l k ) (i = 6, · · · , 14). The relative deviation of intracellular substances with regard to parameter vector p against the disturbance p is defined as For given parameter vector p ∈ D ps (j, l k ), the expectation and variance of the relative deviation RSD i (p , p ; j, l k ), i = 6, · · · 14, caused by p uniformly distributed in field of B(p, δ) are respectively defined by where w 1 , w 2 are the weight coefficients for the balance of the magnitude.
The identification problem of finding the corresponding parameters such that the computational values of extracellular steady states approximate the experimental results as possible as they can and that the system is the most robust, can be formulated as: where Eq.(34) is said to be continuous state inequality constraints [29], which must be satisfied at every point in the time horizon [t 0 , t f ], Eq.(35) is the state vector x(t; x 0 , p j,k , j, l k ) ∈ S x (j, l k ) to be an approximately stable solution, Eq.(36) is a quality constraint for ensuring an acceptable level of the system NHDS(j, l k ) better fits given experimental data.
Theorem 4.5. For given (j, l k ) ∈ I N × S l , assumed that feasible region of IPM is nonempty, there exists an optimal solution of IPM, that is ∃p * j,k ∈D ps (j, l k ) such that J (p * j,k , l k )≤J (p j,k , l k ), ∀p∈D ps (j, l k ).
Proof. By the Properties 2, 3 and definition of the J (p j,k , l k ) , we know that J (p j,k , l k ) is continuous at p j,k ∈ D ps (j, l k ). Since D ps (j, l k ) is compact in R 31 , there exists an optimal solution of IPM.
In the identification problem IPM, the objective function J (p j,k , l k ) is continuous with respect to the parameter vector p j,k ∈ R 31 (j ∈ I N , k ∈ I 9 ) and discrete with respect to the path parameter l i k (i ∈ I 4 , k ∈ I 9 ). We consider 9 possible paths and N experiments, so in IPM we have 31 × 9 × N continuous variables and 4 × 9 × N discrete variables. Therefore, the burden of calculation for the definition of biological robustness will increase significantly as the number of experiments increases. Through observation, we find that system parameters p j,k and path parameters l k are independent of each other in the objective function J (p j,k , l k ), the original problem thence can be decomposed into two sub-problems to optimize the path parameters and system parameters.
Firstly, given l k ∈ S l , the identification subproblem corresponding to system NHDS(j, l k ) denoted by IPM 1 (l k ), can be formulated as: x(t; x 0 , p, j, l k ) ∈ W ad , ∀t∈[t 0 , t f ], p j,k ∈ D ps (j, l k ).
Since l k ∈ S l is given, there is only continuous parameter vector p j,k ∈ D ps (j, l k ) in IPM 1 (l k ). Secondly, the identification subproblem concerning l k ∈ S l , denoted by IPM 2 , can be formulated as: where p * j,k is the optimal parameter vectors of IPM 1 (l k ). Let l * k := argmin{J (p * j,k , l k ) | l k ∈ S l )}. Then, the optimal robustness value and the optimal discrete parameter vectors of IPM 2 are J (p * j,k , l * k ) and l * k .
5. Numerical solutions. The identification problem IPM given by last section can be treated as a optimal parameter selection problem. Similar with the reference [23], any kind of gradient-based optimization methods can be used to solve the problem IPM. However, those gradient-based algorithm are designed to find local optimal solutions. In order to solve the IPM problem globally, we combined a gradient-based and particle swarm optimization method to design a global search algorithm.

Constraint transcription.
When we intend to solve the problem IPM, our main difficulty lies in solving the subproblem IPM 1 (l k ). In nature, given l k ∈ S l , IPM 1 (l k ) is a constrained optimization problem. However, it is difficult to cope with the continuous state inequality constraints. To surmount these difficulties, let g i (t; p j,k , j, l k ) := x * i − x i (t; x 0 , p j,k , j, l k ), g 14+i (t; p j,k , j, l k ) := x i (t; x 0 , p j,k , j, l k ) − x * i , i ∈ I 14 .
The condition x(t; x 0 , p, j, l k ) ∈ W ad is equivalently transcribed into G(p j,k , j, l k ) = 0, where G(p j,k , j, l k ) := 28 i=1 t f t0 min{0, g i (t; p j,k , j, l k )}dt. However, G(p j,k , j, l k ) is non-smooth and the non-differentiable points are at g i (t; p j,k , j, l k ) = 0, i ∈ I 28 . Consequently, standard optimization routines would have difficulties with this type of equality constrains. The following smoothing technique introduced in [15] is to replace min{0, g i (t; p j,k , j, l k )} with ϕ ε,i (t; p j,k , j, l k ), where where ε > 0 is an adjustable parameter controlling the accuracy of the approximation. Then, the equality constraint (40) now can be approximated bỹ where γ > 0 is an adjustable parameter. Hence, with this approximation scheme, the concentration bounds in problem IPM and subproblem IPM 1 (l k ) are approximated by the single canonical constraint (41). This constraint is a standard constraint and can be readily handled using the computational algorithm described later in Section 5.2. Then, IPM and IPM 1 (l k ) can be approximated by the following approximate problems called IPM ε,γ and IPM 1,ε,γ (l k ): Theorem 5.1. Given (j, l k ) ∈ I N × S l , for each ε > 0, there exists a γ(ε) > 0 such that for all γ, 0 < γ < γ(ε), any feasible solution of the model IPM 1,ε,γ (l k ), and IPM ε,γ , is also a feasible solution of the model IPM 1 (l k ), and IPM.

Remark 1.
If the parameter p j,k is not feasible, we can move the parameter towards the feasible region in the direction of the gradients of constraintG ε,γ (p j,k , j, l k ) with respect to parameter p j,k . Then, the gradient of the constraintG ε,γ (p j,k , j, l k ) will be needed and the gradient information is given by the next theorem.
Proof. The proof can be completed using the method of Chapter 3 in [4].

Optimization algorithms.
Particle Swarm Optimization (PSO) is an evolutionary computational method which is based on swarm intelligence. PSO was proposed by Kennedy and Eberhart inspired by the research of the artificial life in [11]. Since then, it has been widely used to solve a broad range of optimization problems. PSO is based on a simple concept, similar to evolutionary optimization methods, PSO is a derivative-free, population-based global search algorithm. Therefore, the computation time is short and it requires few memories. It was originally developed for nonlinear optimization problems with continuous variables and discrete variables [1,3,10,26].
To solve the parameters identification, we design a parallel algorithms as follows, is similar to [20].

Parallel Algorithms 1:
Parallel algorithm of solving identification problem IPM 1,γ,ε (l k ), is as follows: • p j,k := (p 1 j,k , p 2 j,k , · · · , p 31 j,k ) T be the kinetic parameter vector to be determined. • M ρ is the total number of particles in the swarm.
• A sup is a sufficient large number.
• V min and V max are vectors containing the minimum and maximum particle velocities. • θ is the iteration index, Θ max is the minimum and maximum number of iterations. • p w best is the personal best position of the wth particle. • p g best is the best position of whole population. • c 1 and c 2 are the cognitive and social scaling parameters.
• r 1 and r 2 are the random numbers with the uniform distribution on [0, 1].
• is the inertia weight, min and max are the minimum and maximum inertia weights.
Step 1: Allocate n ρ processes for the algorithm. Let the population size be M ρ . Then the number of particles is n w = M ρ /n ρ . Initialized data on the root process as follows: Step 1.1: Read known data such as µ m , K s , K mG , K mP , W ad , P ad ; Step 1.3: Read experimental data as D j , C s0,j , y j = (y 1,j , · · · , y 5,j ) T , j ∈ I N .
Step 2: Broadcast data of root process to all processes.
Step 4: Execute the following procedure on each process, and denote the ID of process as ρ.
Step 4.1: Set θ := 1, according to the uniform distribution, randomly generate the positions of n w particles from the region P ad . Denote the position and velocity of particle by p ρ,w,θ j,k ∈ P ad and ν ρ,w,θ j,k ∈ [V min , V max ], w = 1, 2, · · · , n w .
Step 4.5: Choose the particle p ρ,w,θ j,k , (1 ≤ w ≤ n w ), which has the best fitness value in the wth group, if J(p ρ,w,θ j,k , j, l k ) < J(p w best , j, l k ), then p w best := p ρ,w,θ j,k , if J(p w best , j, l k ) < J(p g best , j, l k ), then p g best := p w best .
Step 4.6: If w < n w , then set w := w + 1, and go to Step 4.3.
Step 5: Gather J(p g best , j, l k ), p g best from process ρ into the root process.
Step 6: Execute the following procedure on each process.
Step 6.1: Set θ := θ + 1. Change the velocity and position of each particle then return Step 4.3, update particles according to the following rules: where ν ρ,w,θ j,k (i) and p ρ,w,θ j,k (i) represent the ith (i ∈ I 31 ) component of velocity and position of the wth (w ∈ I nw ) particle of the ρth (ρ ∈ I nρ ) process at the θ iteration, respectively. Change (θ) according to the following formula: When the parameters are out of their bounds, they are treated as follows: Step 6.2: If θ < Θ max , then set θ := θ + 1 and go to Step 4.2, Otherwise go to Step 7.
Step 7: Execute the following procedure on the root process.
Step 7.1: Output the whole population p g best , and compute J(p g best , j, l k ).
Step 7.2: If j < N , then set j := j + 1, and go to Step 4. Otherwise, compute J (p g best , l k ) := 1 N n j=1 J(p g best , j, l k ), then stop.
Combining Algorithm 1 with Theorem 5.1, we can develop the following algorithm to solve the IPM 1 (l k ).
Step 6: Output p * j,k,ε,γ and l k , then stop. Remark 2. In Algorithm 2, ε is a parameter controlling the accuracy of the smoothing approximation. γ is a parameter controlling the feasibility of the constraint (41). ε andγ are two predefined parameters ensuring the termination of the algorithm.
Remark 3. It is important for the validity of Algorithm 2 to choose the parameters α, β,ε andγ. Especially, the parameters α and β must be chosen less than 1.ε andγ are two sufficiently small values such that Algorithm 2 can be terminated.
Remark 4. The formula of the E i (t s ; RSD i , p, j, l k ) is expressed as an integral form in Eq. (30), which would be difficult to compute directly. Hence, we need divide with the uniform distribution on the perturbation space m 0 subdivisions. The formula of the E i (t s ; RSD i , p, j, l k ) can be change as follows: where p q (q ∈ I m0 ) are uniformly distributed in field of B(p, δ) , p ∈ D ps (j, l k ), and δ is the radius of this field.
To find a appropriate value of the number of subdivisions m 0 , we propose the following algorithm: Algorithm 3: Step 1 : Given τ > 0, choose two positive integers A < B, and set δ := 0.
Step 2 : It follows from (48) that we compute the value of E A (t s ; RSD i , p, j, l k ) and E B (t s ; RSD i , p, j, l k ).
Step 4 : Seek for the appropriate value of the number of subdivisions m 0 is using by Algorithm 3. According to Eqs.(30)- (32) and Eq.(48), the results for m 0 ≥ 5000 show far less variation in the value of the cost function than the results for m 0 < 5000. That is, we can see that 5000 is an appropriate value of m 0 .
The identification model IPM includes 27 nonlinear dynamical systems, 837 continuous variables and 108 discrete variables. Problem IPM is solved by using Algorithms 1 on Lenovo LingYun Large-scale Clusters Supercomputer (which is located in School of Dalian University of Technology, China). As shown in Table 3, for the metabolic systems NHDS(l k ), k = 2, 4, 7, their robustness index are set to be +∞, because the sets D ps (j, l k ), k = 2, 4, 7 are empty for all j ∈ I N even if the number of random samples is large enough, which indicates that these metabolic systems (NHDS(l 2 ), NHDS(l 4 ) and NHDS(l 7 ) ) are all very poor in representing the actual fermentation process, the optimal performance index of IPM with l 9 = (1, 1, 1, 1) T is the smallest, which demonstrates that it is most possible that 1,3-PD and glycerol passes the cell membrane by passive diffusion coupling with active transport. The corresponding optimal dynamical parameter vector p * (l 9 ) with l 9 = (1, 1, 1, 1) T is shown in Table 4. The approximate stable state vector x * under the optimal dynamical parameter vector p * (l 9 ) is shown in Table 5. The comparison of the three extracellular (biomass, glycerol, 1,3-PD) concentrations between experimental data and computational results is shown in Fig. 2, and Fig. 3 -Fig. 4 show the simulated results of intracellular (glycerol, 3-HPA,1,3-PD, mR, R, mGDHt, DGHt, mPDOR, PDOR) concentrations, respectively.

Discussions and conclusions.
In this paper, we establish a complex metabolic system composed of 9 possible pathways and formulated the corresponding nonlinear hybrid dynamical to describe the continuous culture fermentation of glycerol by K. pneumoniae. We then discuss  some basic properties of the metabolic system. A quantitative biological robustness for the concentrations of intracellular substances is defined by penalizing a weighted sum of the expectation and variance of the relative deviation between system outputs before and after parameters are perturbed. In order to determine the parameters, a parameter identification model is presented. The identifiability of the parameter identification model is also proved. Because of the parameter variables of the system concluding continuous and discrete variables, and the parameter identification model does not has analytical solution, we construct a parallel optimization algorithm, based on particle swarm optimization migration (PSO), for solving the identification model. Numerical results indicate that the nonlinear dynamical system is suitable for the process of continuous culture. The results in this work are consistent with previous conjectures in fermentation.
In a future work, our efforts will focus on multistage analysis and optimal control of glycerol input in continuous culture. Additionally, necessary condition of optimality for the system identification model is waiting for publication in another paper.