# Modeling and Parameter Identification Involving 3-Hydroxypropionaldehyde Inhibitory Effects in Glycerol Continuous Fermentation

Authors: Zhaohua Gong; Chongyang Liu; Yongsheng Yu
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Mathematical modeling and parameter estimation are critical steps in the optimization of biotechnological processes. In the 1,3-propanediol (1,3-PD) production by glycerol fermentation process under anaerobic conditions, 3-hydroxypropionaldehyde (3-HPA) accumulation would arouse an irreversible cessation of the fermentation process. Considering 3-HPA inhibitions to cells growth and to activities of enzymes, we propose a novel mathematical model to describe glycerol continuous cultures. Some properties of the above model are discussed. On the basis of the concentrations of extracellular substances, a parameter identification model is established to determine the kinetic parameters in the presented system. Through the penalty function technique combined with an extension of the state space method, an improved genetic algorithm is then constructed to solve the parameter identification model. An illustrative numerical example shows the appropriateness of the proposed model and the validity of optimization algorithm. Since it is difficult to measure the concentrations of intracellular substances, a quantitative robustness analysis method is given to infer whether the model is plausible for the intracellular substances. Numerical results show that the proposed model is of good robustness.


---

## Body

## 1. Introduction

Microbial conversion of glycerol to 1,3-propanediol (1,3-PD) is particularly attractive in that the process is relatively easy and does not generate toxic byproducts. 1,3-PD has numerous applications in polymers, cosmetics, foods, lubricants, and medicines. Industrial 1,3-PD production has attracted attention as an important monomer to synthesize a new type of polyester, polytrimethylene terephthalate (PTT) [1]. However, compared with chemical routes, microbial production is difficult to obtain a high 1,3-PD concentration. The fermentation of glycerol by Klebsiella pneumoniae (K. pneumoniae) under anaerobic conditions is summarized in Figure 1. In the reductive pathway, 3-hydroxypropionaldehyde (3-HPA) is a toxic intermediary metabolite and its accumulation would arouse inhibitions to cells growth and to activities of the enzymes (such as glycerol dehydratase (GDHt) and 1,3-propanediol oxidoreductase (PDOR)) in glycerol metabolism [2–5].Figure 1
Anaerobic metabolism pathways of glycerol byK. pneumoniae [6].It is critical to formulate the fermentation process using a precise mathematical model in the optimization of biotechnological processes. An excess kinetic model for substrate consumption and product formation was established in previous studies [7–10]. The models have been studied for parameter identification [11] and optimal control [12–14] in fed-batch fermentation process. However, the intermediate and intracellular substances or enzymes of glycerol metabolism are not taken into consideration in those models. In fact, some important intermediate substances (such as 3-HPA), intracellular substances (such as 1,3-PD), and enzymes GDHt and PDOR play significant roles in glycerol metabolism. A mathematical model of glycerol fermentation concerning enzyme-catalytic reductive pathway and transports of glycerol and 1,3-PD across cell membrane was established in [6]. Although the achieved results are interesting, the effect of 3-HPA on cells growth is ignored. Moreover, that model is based on an assumption that 3-HPA inhibits the activities of the enzymes GDHt and PDOR all the time. In fact, there exists inhibitory effect of 3-HPA on cells growth owing to its toxicity. In addition, only when the accumulation of 3-HPA reaches some critical concentration, the inhibitions to enzymes can occur [3, 4].Robustness is one of the fundamental characteristics of biological systems. By saying that a system is robust we imply that a particular function or characteristic of the system is preserved despite changes in the operating environment [15]. For robust biological systems, we expect that mathematical models attempting to explain these systems should also be robust [16]. In this paper, we are interested in the robustness to variations in kinetic parameters and use it to validate the plausibility of the mathematical model. This topic has been studied by the sensitivity analysis technique [17–19], that is, repeated simulations by varying one parameter while holding all others fixed. However, single parameter insensitivity may not be sufficient owing to interactions between several parameters. Therefore, new methods are needed for studying multiparameter robustness.Considering 3-HPA inhibitions to cells growth and to activities of the enzymes GDHt and PDOR in glycerol metabolism, we propose a novel mathematical model to describe 1,3-PD production byK. pneumoniae in continuous cultures. Some properties of the model, such as existence and uniqueness of the solution, continuity of the solution in kinetic parameters, and compactness of the set of feasible parameters, are discussed. Furthermore, a parameter identification model is established to determine the kinetic parameters in the presented system. Basing on the penalty function technique and an extension of the state space method, an improved genetic algorithm (GA) is then constructed to solve the identification model. Numerical example shows the appropriateness of the proposed system and the validity of optimization algorithm. Finally, a quantitative robustness analysis method is given to infer whether the model is robust, and numerical result shows that the proposed system is of good robustness.This paper is organized as follows. In Section2, the kinetic model is formulated to describe continuous fermentation process, whose important properties are also discussed. In Section 3, a parameter identification model is presented and an optimization algorithm is developed. Section 4 explores the robustness analysis of the proposed dynamical system. Finally, conclusions are provided.


## 2. Mathematical Model in Continuous Culture and Its Properties

### 2.1. Mathematical Model

During continuous fermentation of glycerol metabolism byK. pneumoniae under anaerobic conditions, glycerol is fed to the reactor continuously. As stated in [6], we assume that(H1)
the transport of extracellular glycerol across cells membrane by passive diffusion and by glycerol transport facilitator;(H2)
intracellular 1,3-PD is expected to be diffused from the intracellular environment to the extracellular medium in the fermentative broth.Letx(t):=(x1(t),x2(t),x3(t),x4(t),x5(t),x6(t),x7(t),x8(t))∈R8 be the state vector, and xi(t), i=1,2…,8, respectively denote the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, extracellular acetic, extracellular ethanol, intracellular glycerol, intracellular 3-HPA, and intracellular 1,3-PD at time t in reactor. D is the dilution rate, and Cs0 is the initial glycerol concentration in feed.Since 3-HPA is a toxic intermediary metabolite, its inhibition to the specific cellular growth rate is introduced besides substrate and products inhibitions. Therefore, the specific cellular growth rateμ [6] is modified as
(2.1)μ=μmx2x2+ks*(1-x2x2*)(1-x3x3*)(1-x4x4*)(1-x5x5*)(1-x7x7*).
Under anaerobic conditions at 37°C and pH 7.0, the maximum specific growth rate, μm, and Monod constant, ks*, are 0.67h-1 and 0.28 mmol L-1, respectively. The critical concentrations xi*, i=2,3,4,5,7, are 2039, 1036, 1026, 360.9, and 300 mmolL-1 [6], respectively.The specific consumption rate of substrateq2, and the specific formation rates of products q4 and q5 are expressed by the following equations based on previous works [6, 8, 10]
(2.2)q2=k6+μk7+k8x2x2+k9,q4=-k10+μk11+k12x2x2+k13,q5=k14+μk15+k16x2x2+k17.The governing equations based on mass balance [6, 10] can still be used to describe the concentrations of substrate and products for glycerol metabolism, for example, biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol, and intracellular glycerol which are described by (2.3)–(2.7), respectively
(2.3)x˙1=(μ-D)x1,(2.4)x˙2=D(Cs0-x2)-q2x1,(2.5)x˙3=k1(x8-x3)x1-Dx3,(2.6)x˙i=qix1-Dxi,i=4,5,(2.7)x˙6=1k2(k3x2x2+k4+k5(x2-x6)-q2)-μx6.The reductive pathway is emphasized because 3-HPA is the key intermediate for 1,3-PD production. 3-HPA accumulation during fermentation process can cause growth cessation and low product formation [2, 4]. Moreover, when the accumulation of 3-HPA reaches some critical concentration, the inhibitions to enzymes GDHt and PDOR can occur [3, 4]. So the intracellular concentration change of 3-HPA can be described by
(2.8)x˙7=k18u1x6km1*(1+IR+(x7-a*)((x7-a*)/k19))+x6-k20u2x7km2*+x7(1+IR+(x7-a*)((x7-a*)/k21))-μx7,
where km1*, km2* are 0.53 and 0.14 mmol L-1 [20, 21], respectively. a* is the critical concentration of 3-HPA beyond which the inhibitions to the activities of GDHt and PDOR occur. Moreover, in (2.8),
(2.9)IR+(z)={1,z>0,0,z≤0.u1 and u2 are the specific activities of GDHt and PDOR in vitro, which can be described by the following equations:
(2.10)u1=k22-k23μ-k24x7x7+k25,u2=k26-k27μ-k28x7x7+k29.The intracellular 1,3-PD concentration depends on the conversion of 3-HPA catalyzed by PDOR whose activity is inhibited by the substrate, the diffusion from the intercellular to the extracellular and the dilution effect on cell growth, so whose variation can be formulated by(2.11)x˙8=k20u2x7km2*+x7(1+IR+(x7-a*)((x7-a*)/k21))-k30(x8-x3)-μx8.Now, letk:=(k1,k2,…,k30) be the kinetic parameter vector to be identified. Denote u:=(CS0,D) and let the right-hand sides of (2.3)–(2.7), (2.8), and (2.11) be f(x,u,k). Then, the proposed mathematical model can be rewritten as the following nonlinear dynamical system:
(2.12)x˙=f(x,u,k)x(0)=x0,t∈[0,T],
where T is the steady-state moment of the continuous fermentation process.


### 2.2. Properties of the Dynamical System

To begin with, we introduce some symbols which will be used below. LetK⊂R30 be the admissible set of the kinetic parameter vector k. Let x*:=(0.001,100,0,0,0,0,0,0) and x*:=(5,2039,1036,1026,360.9,300,2039,1036) denote the lower and upper bounds of the state vector x, respectively. Let Wa:=∏i=18[xi*,xi*]⊂R+8 be the admissible set of x, and let u∈U:=[110.96,1883]×[0.1,0.5] be the admissible set of initial glycerol concentration in feed medium Cs0 and dilution rate D.For the system (2.12), we assume that(H3)
the setK⊂R30 is a nonempty bounded closed set;(H4)
the absolute difference between extracellular and intracellular 1,3-PD and that of glycerol concentration is bounded, that is,∃M1>0 and M2>0 such that
(2.13)|x8(t)-x3(t)|≤M1,∀t∈[0,T],|x2(t)-x6(t)|≤M2,∀t∈[0,T].Under the assumptions (H3) and (H4), we can easily verify the following properties of the velocity vector fieldf(x,u,k).Property 1.
For anyk∈K and u∈U, the function f(x,u,k) is locally Lipschitz continuous in x on Wa.Property 2.
For anyk∈K and u∈U, the function f(x,u,k) satisfies linear growth condition, that is, there exist constants α,β>0 such that
(2.14)∥f(x,u,k)∥≤α+β∥x∥,∀x∈Wa,
where ∥·∥ is Euclidean norm.Proof.
For givenk∈K and u∈U, let C1:=μm, C2:=|k6|+μm/|k7|+|k8|+D, C3:=|k1|M1, C4:=|k10|+μm|k11|+|k12|, C5:=|k14|+μm|k15|+|k16|, C6:=(1/|k2|)(|k3|+|k5|M2+C2-D), C7:=|k18|(|k22|+|k23|μm+|k24|)+|k20|(|k26|+|k27|μm+|k28|), C8:=|k20|(|k26|+|k27|μm+|k28|+|k30|M1. Let Li:=Ci+D, i=1,3,4,5, L2:=max{DCs0,C2}, Li:=max{Ci,μm}, i=6,7,8. Then, we can obtain
(2.15)|f1(x,u,k)|≤(μ+D)|x1|≤L1(∥x∥+1),|f2(x,u,k)|≤DCs0+D|x2|+|q2||x1|≤L2(∥x∥+1),|f3(x,u,k)|≤M2|k1x1|+D|x3|≤L3(∥x∥+1),|fi(x,u,k)|≤|qi||x1|+D|xi|≤Li(∥x∥+1),i=4,5,|fi(x,u,k)|≤Ci+μxi≤Li(∥x∥+1),i=6,7,8.
Finally, set L′:=max{Li,i∈I8:={1,2,…,8}}, then we have (2.14) holds with α=β=22L′. The proof is completed.Then, the existence and uniqueness of the solution for the system (2.12) can be confirmed in the following theorem.Theorem 2.1.
For anyk∈K and u∈U, the system (2.12) with given initial state x0∈Wa has a unique solution denoted by x(·;u,k). Moreover, x(t;u,k) is continuous in k on K.
(2.16)x(t;u,k)=x0+∫0tf(x(s),u,k)ds,∀t∈[0,T].Proof.
The proof can be obtained from Properties1, and 2 and the theory of ordinary differential equations [22].Givenx0∈Wa, we define the solution set 𝒮0 of the system (2.12) as follows:
(2.17)𝒮0≔{x(·;u,k)∣x(t;u,k)isasolutionofthesystem(2.12)correspondingtok∈Kandu∈Uforanyt∈[0,T]}.
Since the concentrations of biomass, glycerol, 3-HPA, and products are restricted in Wa during the actual continuous cultures, the set of admissible solutions is
(2.18)𝒮:={x(·;u,k)∈𝒮0∣x(t;u,k)∈Wa∀t∈[0,T]}.
Furthermore, let the set of feasible parameter vectors corresponding to 𝒮 be
(2.19)ℱ:={k∈K∣x(·;u,k)∈𝒮}.From Theorem2.1 and the above definitions, we have the following result.Theorem 2.2.
The feasible parameter setℱ defined by (2.19) is a compact set.Proof.
In view of the compactness ofK, we obtain that ℱ is a bounded set. Moreover, for any sequence {ki}i=1∞∈K, there exists at least a subsequence {k^i}⊆{ki} such that k^i→k^ as i→∞. It follows from Theorem 2.1 that x(·;u,k^i)∈𝒮0 and x(·;u,k^)∈𝒮0. Since Wa is a compact set, we must conclude that x(·;u,k^)∈𝒮, which implies the closeness of set ℱ. The proof is completed.


## 2.1. Mathematical Model

During continuous fermentation of glycerol metabolism byK. pneumoniae under anaerobic conditions, glycerol is fed to the reactor continuously. As stated in [6], we assume that(H1)
the transport of extracellular glycerol across cells membrane by passive diffusion and by glycerol transport facilitator;(H2)
intracellular 1,3-PD is expected to be diffused from the intracellular environment to the extracellular medium in the fermentative broth.Letx(t):=(x1(t),x2(t),x3(t),x4(t),x5(t),x6(t),x7(t),x8(t))∈R8 be the state vector, and xi(t), i=1,2…,8, respectively denote the concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, extracellular acetic, extracellular ethanol, intracellular glycerol, intracellular 3-HPA, and intracellular 1,3-PD at time t in reactor. D is the dilution rate, and Cs0 is the initial glycerol concentration in feed.Since 3-HPA is a toxic intermediary metabolite, its inhibition to the specific cellular growth rate is introduced besides substrate and products inhibitions. Therefore, the specific cellular growth rateμ [6] is modified as
(2.1)μ=μmx2x2+ks*(1-x2x2*)(1-x3x3*)(1-x4x4*)(1-x5x5*)(1-x7x7*).
Under anaerobic conditions at 37°C and pH 7.0, the maximum specific growth rate, μm, and Monod constant, ks*, are 0.67h-1 and 0.28 mmol L-1, respectively. The critical concentrations xi*, i=2,3,4,5,7, are 2039, 1036, 1026, 360.9, and 300 mmolL-1 [6], respectively.The specific consumption rate of substrateq2, and the specific formation rates of products q4 and q5 are expressed by the following equations based on previous works [6, 8, 10]
(2.2)q2=k6+μk7+k8x2x2+k9,q4=-k10+μk11+k12x2x2+k13,q5=k14+μk15+k16x2x2+k17.The governing equations based on mass balance [6, 10] can still be used to describe the concentrations of substrate and products for glycerol metabolism, for example, biomass, extracellular glycerol, extracellular 1,3-PD, acetate, ethanol, and intracellular glycerol which are described by (2.3)–(2.7), respectively
(2.3)x˙1=(μ-D)x1,(2.4)x˙2=D(Cs0-x2)-q2x1,(2.5)x˙3=k1(x8-x3)x1-Dx3,(2.6)x˙i=qix1-Dxi,i=4,5,(2.7)x˙6=1k2(k3x2x2+k4+k5(x2-x6)-q2)-μx6.The reductive pathway is emphasized because 3-HPA is the key intermediate for 1,3-PD production. 3-HPA accumulation during fermentation process can cause growth cessation and low product formation [2, 4]. Moreover, when the accumulation of 3-HPA reaches some critical concentration, the inhibitions to enzymes GDHt and PDOR can occur [3, 4]. So the intracellular concentration change of 3-HPA can be described by
(2.8)x˙7=k18u1x6km1*(1+IR+(x7-a*)((x7-a*)/k19))+x6-k20u2x7km2*+x7(1+IR+(x7-a*)((x7-a*)/k21))-μx7,
where km1*, km2* are 0.53 and 0.14 mmol L-1 [20, 21], respectively. a* is the critical concentration of 3-HPA beyond which the inhibitions to the activities of GDHt and PDOR occur. Moreover, in (2.8),
(2.9)IR+(z)={1,z>0,0,z≤0.u1 and u2 are the specific activities of GDHt and PDOR in vitro, which can be described by the following equations:
(2.10)u1=k22-k23μ-k24x7x7+k25,u2=k26-k27μ-k28x7x7+k29.The intracellular 1,3-PD concentration depends on the conversion of 3-HPA catalyzed by PDOR whose activity is inhibited by the substrate, the diffusion from the intercellular to the extracellular and the dilution effect on cell growth, so whose variation can be formulated by(2.11)x˙8=k20u2x7km2*+x7(1+IR+(x7-a*)((x7-a*)/k21))-k30(x8-x3)-μx8.Now, letk:=(k1,k2,…,k30) be the kinetic parameter vector to be identified. Denote u:=(CS0,D) and let the right-hand sides of (2.3)–(2.7), (2.8), and (2.11) be f(x,u,k). Then, the proposed mathematical model can be rewritten as the following nonlinear dynamical system:
(2.12)x˙=f(x,u,k)x(0)=x0,t∈[0,T],
where T is the steady-state moment of the continuous fermentation process.


## 2.2. Properties of the Dynamical System

To begin with, we introduce some symbols which will be used below. LetK⊂R30 be the admissible set of the kinetic parameter vector k. Let x*:=(0.001,100,0,0,0,0,0,0) and x*:=(5,2039,1036,1026,360.9,300,2039,1036) denote the lower and upper bounds of the state vector x, respectively. Let Wa:=∏i=18[xi*,xi*]⊂R+8 be the admissible set of x, and let u∈U:=[110.96,1883]×[0.1,0.5] be the admissible set of initial glycerol concentration in feed medium Cs0 and dilution rate D.For the system (2.12), we assume that(H3)
the setK⊂R30 is a nonempty bounded closed set;(H4)
the absolute difference between extracellular and intracellular 1,3-PD and that of glycerol concentration is bounded, that is,∃M1>0 and M2>0 such that
(2.13)|x8(t)-x3(t)|≤M1,∀t∈[0,T],|x2(t)-x6(t)|≤M2,∀t∈[0,T].Under the assumptions (H3) and (H4), we can easily verify the following properties of the velocity vector fieldf(x,u,k).Property 1.
For anyk∈K and u∈U, the function f(x,u,k) is locally Lipschitz continuous in x on Wa.Property 2.
For anyk∈K and u∈U, the function f(x,u,k) satisfies linear growth condition, that is, there exist constants α,β>0 such that
(2.14)∥f(x,u,k)∥≤α+β∥x∥,∀x∈Wa,
where ∥·∥ is Euclidean norm.Proof.
For givenk∈K and u∈U, let C1:=μm, C2:=|k6|+μm/|k7|+|k8|+D, C3:=|k1|M1, C4:=|k10|+μm|k11|+|k12|, C5:=|k14|+μm|k15|+|k16|, C6:=(1/|k2|)(|k3|+|k5|M2+C2-D), C7:=|k18|(|k22|+|k23|μm+|k24|)+|k20|(|k26|+|k27|μm+|k28|), C8:=|k20|(|k26|+|k27|μm+|k28|+|k30|M1. Let Li:=Ci+D, i=1,3,4,5, L2:=max{DCs0,C2}, Li:=max{Ci,μm}, i=6,7,8. Then, we can obtain
(2.15)|f1(x,u,k)|≤(μ+D)|x1|≤L1(∥x∥+1),|f2(x,u,k)|≤DCs0+D|x2|+|q2||x1|≤L2(∥x∥+1),|f3(x,u,k)|≤M2|k1x1|+D|x3|≤L3(∥x∥+1),|fi(x,u,k)|≤|qi||x1|+D|xi|≤Li(∥x∥+1),i=4,5,|fi(x,u,k)|≤Ci+μxi≤Li(∥x∥+1),i=6,7,8.
Finally, set L′:=max{Li,i∈I8:={1,2,…,8}}, then we have (2.14) holds with α=β=22L′. The proof is completed.Then, the existence and uniqueness of the solution for the system (2.12) can be confirmed in the following theorem.Theorem 2.1.
For anyk∈K and u∈U, the system (2.12) with given initial state x0∈Wa has a unique solution denoted by x(·;u,k). Moreover, x(t;u,k) is continuous in k on K.
(2.16)x(t;u,k)=x0+∫0tf(x(s),u,k)ds,∀t∈[0,T].Proof.
The proof can be obtained from Properties1, and 2 and the theory of ordinary differential equations [22].Givenx0∈Wa, we define the solution set 𝒮0 of the system (2.12) as follows:
(2.17)𝒮0≔{x(·;u,k)∣x(t;u,k)isasolutionofthesystem(2.12)correspondingtok∈Kandu∈Uforanyt∈[0,T]}.
Since the concentrations of biomass, glycerol, 3-HPA, and products are restricted in Wa during the actual continuous cultures, the set of admissible solutions is
(2.18)𝒮:={x(·;u,k)∈𝒮0∣x(t;u,k)∈Wa∀t∈[0,T]}.
Furthermore, let the set of feasible parameter vectors corresponding to 𝒮 be
(2.19)ℱ:={k∈K∣x(·;u,k)∈𝒮}.From Theorem2.1 and the above definitions, we have the following result.Theorem 2.2.
The feasible parameter setℱ defined by (2.19) is a compact set.Proof.
In view of the compactness ofK, we obtain that ℱ is a bounded set. Moreover, for any sequence {ki}i=1∞∈K, there exists at least a subsequence {k^i}⊆{ki} such that k^i→k^ as i→∞. It follows from Theorem 2.1 that x(·;u,k^i)∈𝒮0 and x(·;u,k^)∈𝒮0. Since Wa is a compact set, we must conclude that x(·;u,k^)∈𝒮, which implies the closeness of set ℱ. The proof is completed.


## 3. Parameter Identification Model and Optimization Algorithm

### 3.1. Parameter Identification Model

Now, we determine the kinetic parameterk in K by constructing an identification problem as follows.Letl be the total number of experiments carried out under different dilution rates and initial glycerol concentrations. For given uj=(Cs0j,Dj), j∈Il, we have the experimental steady-state data of extracellular substances in continuous cultures. Denote the steady-state concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, and ethanol as y1j,y2j,y3j,y4j,y5j, correspondingly. Let yj:=(y1j,y2j,y3j,y4j,y5j)∈R+5,j∈Il.In particular, the stable state of a dynamical system (2.12) in the following identification problem is actually referred to the approximate stability defined as follows.Definition 3.1.
For givenuj∈U, a state vector x(Tj;uj,k) is said to be an approximately stable solution within a precision ε>0 if there exists k∈K such that x(·;uj,k) is a solution of the system (2.12) satisfying
(3.1)∥f(x(Tj),uj,k)∥<ε,
where ∥·∥ is the Euclidean norm, and Tj=inf{tε:∥f(x(t),uj,k)∥<εfor allt∈[tε,T]}.Since the orders of magnitude for concentrations involved are different, we adopt the average relative error between the computational valuesx(Tj;uj,k) and the experimental data yj at steady-state moments Tj,j∈Il, as the criterion
(3.2)J(k):=15∑m=15∑j=1l|xm(Tj;uj,k)-ymj|∑j=1lymj.To determine the parameter values of the system (2.12), a parameter identification problem, in which J(k) is taken as the cost function, can be formulated as
(PIP)minJ(k)s.t.∥f(x(Tj),uj,k)∥<ε,j∈Il,k∈ℱ.Following the above properties, we can conclude the following theorem.Theorem 3.2.
Foruj∈U, j∈Il, there exists k*∈ℱ in (PIP) such that J(k*)≤J(k),for allk∈ℱ.


### 3.2. Optimization Algorithm

In(PIP), the constraint k∈ℱ actually involves the constraint of continuous state, that is, x*≤x(t;uj,k)≤x*,t∈[0,Tj],k∈K,j∈Il. In this section, we will develop a computational method for solving our proposed parameter identification problem (PIP). By means of the penalty function technique and an extension of state space method, we transcribe (PIP) into an optimization problem only with box constraint. First of all, we introduce a new state variable z satisfying
(3.3)z˙=g(x(t)),t∈[0,Tj],j∈Ilz(0)=0,
where g(x(t)):=[x(t)-x*]++[x*-x(t)]+, and [x~]+:=∑i=18(max{x~i,0})2.Obviously, the processx(t;uj,k) satisfies the constraint of continuous state if and only if z(Tj)=0,j∈Il. Furthermore, denote the cost function in (PIP)  by
(3.4)J~λ,μ(k):=J(k)+λ∑j=1l∥f(x(Tj),uj,k)∥+μ∑j=1lz(Tj),
where λ and μ are penalty factors. Then (PIP)  can be rewritten as
(PIP(λ,μ))minJ~λ,μ(k)s.t.k∈K.As a result,(PIP(λ,μ)) is an optimization problem only with the box constraint and equivalent to (PIP) as λ→+∞ and μ→+∞.SinceJ~λ,μ(k) is nondifferentiable, we construct an improved genetic algorithm (GA) to solve (PIP(λ,μ)) taking advantage of the problem’s characteristic. Let k~:=(k1,k2,…,k30)∈K. In the improved genetic algorithm, we take k~ as the individual, and (3.4) as the fitness function. Now, we describe the algorithm in detail as follows.Algorithm 3.3.
We have the following steps:Step1. Initialize population size N, the maximal iterations H, penalty factors λ and μ, precision ε, parameters N1,α1, and α2. Set h=0 and randomly generate N individuals by uniform distribution, that is, k~i(h)∈K,i=1,2,…N. Let P(h)={k~i(h)∣i=1,2,…,N}.Step2. Compute the fitness value J~λ,μ(k~i(h)) by taking k~i(h) into the systems (2.12) and (3.3).Step3. Generate the crossover offspring dm(h)(m=1,2,…,N1) by arithmetic crossover and evaluate their fitness values J~λ,μ(dm(h)).Step4. The individual produced by crossover is operated with normality variation until the produced one is in K. The mutation offspring is denoted by gm(h)(m=1,2,…,N1). Compute their fitness values J~λ,μ(gm(h)).Step5. Form the next generation P(h+1) by selecting the best N individuals from the N+2N1 ones.Step6. Set h=h+1,λ=α1λ, and μ=α2μ, if either h>H or no progress is made in the last generations, then output the best individual and stop, otherwise go to Step 3.


### 3.3. Numerical Results

According to the actual continuous fermentation process, the initial statex0=(0.1115gL-1, 495mmolL-1,0,0,0,0,0,0) and 22 groups of the experimental steady-state data are used. Here, 10 groups of experimental steady-state data under substrate-limited conditions and 12 groups of experimental steady-state data under substrate-sufficient conditions. The critical value a* is taken as the value in [23]. Moreover, the admissible set of parameter vectors K is taken as the decrements and increments of 0.9 times the kinetic parameter values in [6]. By applying Algorithm 3.3, we obtained the optimal parameter vector k* of the system (2.12) under the substrate-limited and the substrate-sufficient conditions shown in Table 1. Accordingly, the cost function values, respectively, are 1.204×10-2 and 1.585×10-2 under the above two cases. Here, all the computations are performed in Visual C++ 6.0 and numerical results are plotted by MATLAB 7.10.0. In particular, the ODEs in the computation process are numerically calculated by improved Euler method [24] with the relative error tolerance 10-4. The parameters used in Algorithm 3.3 are H=1000, N=50, N1=30, λ=μ=1, ε=0.001, and α1=α2=1.1, respectively. It should be noted that these parameters are derived empirically after numerous experiments. The comparison of three extracellular substances concentrations, that is, biomass, extracellular glycerol, and extracellular 1,3-PD, between experimental steady-state data and computational results under substrate-limited conditions are shown in Figures 2, 3, and 4. Furthermore, The comparisons of three extracellular substance concentrations, that is, biomass, extracellular glycerol, and extracellular 1,3-PD, between experimental steady-state data and computational results under substrate-sufficient conditions are also shown in Figures 5, 6, and 7. From the above figures, we can see that the simulation results can approximate the experimental steady-state data well. Thus, the mathematical model considering 3-HPA inhibitions to cells growth and to activities of enzymes can well describe the continuous fermentation process.Table 1
The optimal parameter vector in the system (2.12).
Conditions
Optimal parameter vectork∗
(16.93, 0.015, 375.9, 0.134, 9.5, 162.8,
0.013, 0.014, 0.518, 0.135, 6.239, 203.6,
Substrate-limited condition
37.5, 0.042, 18.43, 3.13, 0.006489,18.83,
12.82, 1.843, 52.03, 10.91, 5.045, 2.983,
0.567, 0.087, 1.852, 6.225, 0.346, 0.138)
(36.37, 0.287, 456, 0.878, 0.5, 162.8,
0.036, 0.087, 2.235, 2.216, 11.25, 282.1,
Substrate-sufficient condition
39.12, 0.629, 9.841, 3.854, 0.006678,
13.09, 10.77, 0.590, 62.83, 7.833, 17.39,
5.624, 0.353, 0.6145, 0.836, 8.094, 1.254, 0.138)Figure 2
The concentrations of biomass change with respect to fermentation time under substrate-limited conditions foru=(295,0.25).Figure 3
The concentrations of glycerol change with respect to fermentation time under substrate-limited conditions foru=(295,0.25).Figure 4
The concentrations of 1,3-PD change with respect to fermentation time under substrate-limited conditions foru=(295,0.25).Figure 5
The concentrations of biomass change with respect to fermentation time under substrate-sufficient conditions foru=(330,0.4).Figure 6
The concentrations of glycerol change with respect to fermentation time under substrate-sufficient conditions foru=(330,0.4).Figure 7
The concentrations of 1,3-PD change with respect to fermentation time under substrate-sufficient conditions foru=(330,0.4).


## 3.1. Parameter Identification Model

Now, we determine the kinetic parameterk in K by constructing an identification problem as follows.Letl be the total number of experiments carried out under different dilution rates and initial glycerol concentrations. For given uj=(Cs0j,Dj), j∈Il, we have the experimental steady-state data of extracellular substances in continuous cultures. Denote the steady-state concentrations of biomass, extracellular glycerol, extracellular 1,3-PD, acetate, and ethanol as y1j,y2j,y3j,y4j,y5j, correspondingly. Let yj:=(y1j,y2j,y3j,y4j,y5j)∈R+5,j∈Il.In particular, the stable state of a dynamical system (2.12) in the following identification problem is actually referred to the approximate stability defined as follows.Definition 3.1.
For givenuj∈U, a state vector x(Tj;uj,k) is said to be an approximately stable solution within a precision ε>0 if there exists k∈K such that x(·;uj,k) is a solution of the system (2.12) satisfying
(3.1)∥f(x(Tj),uj,k)∥<ε,
where ∥·∥ is the Euclidean norm, and Tj=inf{tε:∥f(x(t),uj,k)∥<εfor allt∈[tε,T]}.Since the orders of magnitude for concentrations involved are different, we adopt the average relative error between the computational valuesx(Tj;uj,k) and the experimental data yj at steady-state moments Tj,j∈Il, as the criterion
(3.2)J(k):=15∑m=15∑j=1l|xm(Tj;uj,k)-ymj|∑j=1lymj.To determine the parameter values of the system (2.12), a parameter identification problem, in which J(k) is taken as the cost function, can be formulated as
(PIP)minJ(k)s.t.∥f(x(Tj),uj,k)∥<ε,j∈Il,k∈ℱ.Following the above properties, we can conclude the following theorem.Theorem 3.2.
Foruj∈U, j∈Il, there exists k*∈ℱ in (PIP) such that J(k*)≤J(k),for allk∈ℱ.


## 3.2. Optimization Algorithm

In(PIP), the constraint k∈ℱ actually involves the constraint of continuous state, that is, x*≤x(t;uj,k)≤x*,t∈[0,Tj],k∈K,j∈Il. In this section, we will develop a computational method for solving our proposed parameter identification problem (PIP). By means of the penalty function technique and an extension of state space method, we transcribe (PIP) into an optimization problem only with box constraint. First of all, we introduce a new state variable z satisfying
(3.3)z˙=g(x(t)),t∈[0,Tj],j∈Ilz(0)=0,
where g(x(t)):=[x(t)-x*]++[x*-x(t)]+, and [x~]+:=∑i=18(max{x~i,0})2.Obviously, the processx(t;uj,k) satisfies the constraint of continuous state if and only if z(Tj)=0,j∈Il. Furthermore, denote the cost function in (PIP)  by
(3.4)J~λ,μ(k):=J(k)+λ∑j=1l∥f(x(Tj),uj,k)∥+μ∑j=1lz(Tj),
where λ and μ are penalty factors. Then (PIP)  can be rewritten as
(PIP(λ,μ))minJ~λ,μ(k)s.t.k∈K.As a result,(PIP(λ,μ)) is an optimization problem only with the box constraint and equivalent to (PIP) as λ→+∞ and μ→+∞.SinceJ~λ,μ(k) is nondifferentiable, we construct an improved genetic algorithm (GA) to solve (PIP(λ,μ)) taking advantage of the problem’s characteristic. Let k~:=(k1,k2,…,k30)∈K. In the improved genetic algorithm, we take k~ as the individual, and (3.4) as the fitness function. Now, we describe the algorithm in detail as follows.Algorithm 3.3.
We have the following steps:Step1. Initialize population size N, the maximal iterations H, penalty factors λ and μ, precision ε, parameters N1,α1, and α2. Set h=0 and randomly generate N individuals by uniform distribution, that is, k~i(h)∈K,i=1,2,…N. Let P(h)={k~i(h)∣i=1,2,…,N}.Step2. Compute the fitness value J~λ,μ(k~i(h)) by taking k~i(h) into the systems (2.12) and (3.3).Step3. Generate the crossover offspring dm(h)(m=1,2,…,N1) by arithmetic crossover and evaluate their fitness values J~λ,μ(dm(h)).Step4. The individual produced by crossover is operated with normality variation until the produced one is in K. The mutation offspring is denoted by gm(h)(m=1,2,…,N1). Compute their fitness values J~λ,μ(gm(h)).Step5. Form the next generation P(h+1) by selecting the best N individuals from the N+2N1 ones.Step6. Set h=h+1,λ=α1λ, and μ=α2μ, if either h>H or no progress is made in the last generations, then output the best individual and stop, otherwise go to Step 3.


## 3.3. Numerical Results

According to the actual continuous fermentation process, the initial statex0=(0.1115gL-1, 495mmolL-1,0,0,0,0,0,0) and 22 groups of the experimental steady-state data are used. Here, 10 groups of experimental steady-state data under substrate-limited conditions and 12 groups of experimental steady-state data under substrate-sufficient conditions. The critical value a* is taken as the value in [23]. Moreover, the admissible set of parameter vectors K is taken as the decrements and increments of 0.9 times the kinetic parameter values in [6]. By applying Algorithm 3.3, we obtained the optimal parameter vector k* of the system (2.12) under the substrate-limited and the substrate-sufficient conditions shown in Table 1. Accordingly, the cost function values, respectively, are 1.204×10-2 and 1.585×10-2 under the above two cases. Here, all the computations are performed in Visual C++ 6.0 and numerical results are plotted by MATLAB 7.10.0. In particular, the ODEs in the computation process are numerically calculated by improved Euler method [24] with the relative error tolerance 10-4. The parameters used in Algorithm 3.3 are H=1000, N=50, N1=30, λ=μ=1, ε=0.001, and α1=α2=1.1, respectively. It should be noted that these parameters are derived empirically after numerous experiments. The comparison of three extracellular substances concentrations, that is, biomass, extracellular glycerol, and extracellular 1,3-PD, between experimental steady-state data and computational results under substrate-limited conditions are shown in Figures 2, 3, and 4. Furthermore, The comparisons of three extracellular substance concentrations, that is, biomass, extracellular glycerol, and extracellular 1,3-PD, between experimental steady-state data and computational results under substrate-sufficient conditions are also shown in Figures 5, 6, and 7. From the above figures, we can see that the simulation results can approximate the experimental steady-state data well. Thus, the mathematical model considering 3-HPA inhibitions to cells growth and to activities of enzymes can well describe the continuous fermentation process.Table 1
The optimal parameter vector in the system (2.12).
Conditions
Optimal parameter vectork∗
(16.93, 0.015, 375.9, 0.134, 9.5, 162.8,
0.013, 0.014, 0.518, 0.135, 6.239, 203.6,
Substrate-limited condition
37.5, 0.042, 18.43, 3.13, 0.006489,18.83,
12.82, 1.843, 52.03, 10.91, 5.045, 2.983,
0.567, 0.087, 1.852, 6.225, 0.346, 0.138)
(36.37, 0.287, 456, 0.878, 0.5, 162.8,
0.036, 0.087, 2.235, 2.216, 11.25, 282.1,
Substrate-sufficient condition
39.12, 0.629, 9.841, 3.854, 0.006678,
13.09, 10.77, 0.590, 62.83, 7.833, 17.39,
5.624, 0.353, 0.6145, 0.836, 8.094, 1.254, 0.138)Figure 2
The concentrations of biomass change with respect to fermentation time under substrate-limited conditions foru=(295,0.25).Figure 3
The concentrations of glycerol change with respect to fermentation time under substrate-limited conditions foru=(295,0.25).Figure 4
The concentrations of 1,3-PD change with respect to fermentation time under substrate-limited conditions foru=(295,0.25).Figure 5
The concentrations of biomass change with respect to fermentation time under substrate-sufficient conditions foru=(330,0.4).Figure 6
The concentrations of glycerol change with respect to fermentation time under substrate-sufficient conditions foru=(330,0.4).Figure 7
The concentrations of 1,3-PD change with respect to fermentation time under substrate-sufficient conditions foru=(330,0.4).


## 4. Robustness Analysis of the Model

For robust biological systems, we expect that mathematical models that attempt to explain these systems should also be robust. Robustness of the model is analyzed in this section. In view of the glycerol dissimilation mechanism, we assume that(H5)
for eachx0∈Wa, k∈ℱ and given u∈U, there exists an approximate stable solution of the system (2.12).


### 4.1. Mathematical Measurement of Robustness

Robustness can be defined as a system’s characteristic that maintains one or more of its functions under external and internal perturbations [15]. In this study, the substance concentrations at steady-state moments are viewed as the quantitative descriptors of the system and the perturbations are the parameters variations in ℱ.LetU- be a finite set whose elements are drawn from U by random distribution. Since the state vector of the system (2.12) takes u∈U as the input parameter, we define the representative steady-state vector as
(4.1)x-(tk;k):=1|U-|∑u∈U-x(tu,k;u,k),
where
(4.2)tu,k:=inf{tε∣∥f(x(t),u,k)∥<εfork∈ℱ,u∈U-,∀t∈[tε,T]},tk:=maxu∈U-{tu,k∣x(tu,k;u,k)isanapproximatelysteady-state vectorof(2.12)withk∈ℱ},
and |U-|denotes  the  cardinal  number  of  the  setU-.LetPS(k,δ):={k~|∥k~-k∥≤δ,k∈ℱ} be the feasible space of the perturbation parameter vector, where δ>0 is sufficiently small positive number to measure the magnitude of the parameter disturbances. Furthermore, let 𝒫 be a perturbation set which is composed of elements randomly generated from PS(k,δ). Denote the set of the representative steady-state vector corresponding to the varied parameters by
(4.3)𝒮-𝒫:={x-(tk;k)∣k∈𝒫}.
For the representative steady-state vectors in 𝒮-𝒫, we are interested in their deviations from the representative steady-state vector corresponding to the optimal parameter vector k*. Define the expectation of the ith component of these deviations as follows:
(4.4)U𝒫i:=1|𝒫|∑k∈𝒫‍(x-i(tk;k)x-i(tk*;k*)-1)2,i∈I8.Based on the above analyses, a mathematical definition of biological robustness can be stated as follows.Definition 4.1.
The robustness measurement of a modelℳ with regard to steady states against a set of parameter perturbations 𝒫 is
(4.5)ℛ𝒫ℳ:=1|Irob|∑i∈Irob1U𝒫i,
where Irob is the state index set involved in robustness analysis.


### 4.2. Algorithm and Numerical Results

On the basis of the definition of robustness, we develop an Algorithm4.2 to compute the robustness to variations of the optimal kinetic parameter vector k*.Algorithm 4.2.
We have the following steps:Step1. Generate N2 perturbation parameter vectors ks∈PS(k*,δ),s=1,2,…,N2 by random distributions domain.Step2. Compute the representative steady-state vectors corresponding to optimal parameters and to perturbation parameters x-(tk*;k*) and x-(tks;ks) according to (4.1).Step3. Evaluate the robustness by the mathematical measurement ℛ𝒫ℳ defined in (4.5).Step4. Repeat Step 1, Step 2, and Step 3 for N3 times and compute the expectation value of the generated sequence {ℛ𝒫ℳ}.According to Algorithm4.2, we have investigated the robustness of the system (2.12) for the intracellular substances. Here, the number of perturbation parameters N2, precision δ, the set Irob, and the repeating times N3 take values 200, 0.05, {6,7,8}, and 100, respectively. Computational results for the substrate-limited and substrate-sufficient cases are listed in Table 2. Furthermore, the simulation curves for intracellular glycerol, 3-HPA and intracellular 1,3-PD under the substrate-limited and the substrate-sufficient conditions are illustrated in Figures 8, 9, 10, 11, 12, and 13. From Table 2 and the numerical results, we can see that the proposed model is of good robustness.Table 2
The robustness of the system (2.12) under substrate-limited and substrate-sufficient conditions.
Conditions
Substrate-limited
Substrate-sufficient
Values
3.97268
2.70685Figure 8
Robustness simulations of intracellular glycerol concentration under substrate-limited conditions foru=(295,0.25).Figure 9
Robustness simulations of 3-HPA concentration under substrate-limited conditions foru=(295,0.25).Figure 10
Robustness simulations of intracellular 1,3-PD concentration under substrate-limited conditions foru=(295,0.25).Figure 11
Robustness simulations of intracellular glycerol concentration under substrate-sufficient conditions foru=(330,0.4).Figure 12
Robustness simulations of 3-HPA concentration under substrate-sufficient conditions foru=(330,0.4).Figure 13
Robustness simulations of intracellular 1,3-PD concentration under substrate-sufficient conditions foru=(330,0.4).


## 4.1. Mathematical Measurement of Robustness

Robustness can be defined as a system’s characteristic that maintains one or more of its functions under external and internal perturbations [15]. In this study, the substance concentrations at steady-state moments are viewed as the quantitative descriptors of the system and the perturbations are the parameters variations in ℱ.LetU- be a finite set whose elements are drawn from U by random distribution. Since the state vector of the system (2.12) takes u∈U as the input parameter, we define the representative steady-state vector as
(4.1)x-(tk;k):=1|U-|∑u∈U-x(tu,k;u,k),
where
(4.2)tu,k:=inf{tε∣∥f(x(t),u,k)∥<εfork∈ℱ,u∈U-,∀t∈[tε,T]},tk:=maxu∈U-{tu,k∣x(tu,k;u,k)isanapproximatelysteady-state vectorof(2.12)withk∈ℱ},
and |U-|denotes  the  cardinal  number  of  the  setU-.LetPS(k,δ):={k~|∥k~-k∥≤δ,k∈ℱ} be the feasible space of the perturbation parameter vector, where δ>0 is sufficiently small positive number to measure the magnitude of the parameter disturbances. Furthermore, let 𝒫 be a perturbation set which is composed of elements randomly generated from PS(k,δ). Denote the set of the representative steady-state vector corresponding to the varied parameters by
(4.3)𝒮-𝒫:={x-(tk;k)∣k∈𝒫}.
For the representative steady-state vectors in 𝒮-𝒫, we are interested in their deviations from the representative steady-state vector corresponding to the optimal parameter vector k*. Define the expectation of the ith component of these deviations as follows:
(4.4)U𝒫i:=1|𝒫|∑k∈𝒫‍(x-i(tk;k)x-i(tk*;k*)-1)2,i∈I8.Based on the above analyses, a mathematical definition of biological robustness can be stated as follows.Definition 4.1.
The robustness measurement of a modelℳ with regard to steady states against a set of parameter perturbations 𝒫 is
(4.5)ℛ𝒫ℳ:=1|Irob|∑i∈Irob1U𝒫i,
where Irob is the state index set involved in robustness analysis.


## 4.2. Algorithm and Numerical Results

On the basis of the definition of robustness, we develop an Algorithm4.2 to compute the robustness to variations of the optimal kinetic parameter vector k*.Algorithm 4.2.
We have the following steps:Step1. Generate N2 perturbation parameter vectors ks∈PS(k*,δ),s=1,2,…,N2 by random distributions domain.Step2. Compute the representative steady-state vectors corresponding to optimal parameters and to perturbation parameters x-(tk*;k*) and x-(tks;ks) according to (4.1).Step3. Evaluate the robustness by the mathematical measurement ℛ𝒫ℳ defined in (4.5).Step4. Repeat Step 1, Step 2, and Step 3 for N3 times and compute the expectation value of the generated sequence {ℛ𝒫ℳ}.According to Algorithm4.2, we have investigated the robustness of the system (2.12) for the intracellular substances. Here, the number of perturbation parameters N2, precision δ, the set Irob, and the repeating times N3 take values 200, 0.05, {6,7,8}, and 100, respectively. Computational results for the substrate-limited and substrate-sufficient cases are listed in Table 2. Furthermore, the simulation curves for intracellular glycerol, 3-HPA and intracellular 1,3-PD under the substrate-limited and the substrate-sufficient conditions are illustrated in Figures 8, 9, 10, 11, 12, and 13. From Table 2 and the numerical results, we can see that the proposed model is of good robustness.Table 2
The robustness of the system (2.12) under substrate-limited and substrate-sufficient conditions.
Conditions
Substrate-limited
Substrate-sufficient
Values
3.97268
2.70685Figure 8
Robustness simulations of intracellular glycerol concentration under substrate-limited conditions foru=(295,0.25).Figure 9
Robustness simulations of 3-HPA concentration under substrate-limited conditions foru=(295,0.25).Figure 10
Robustness simulations of intracellular 1,3-PD concentration under substrate-limited conditions foru=(295,0.25).Figure 11
Robustness simulations of intracellular glycerol concentration under substrate-sufficient conditions foru=(330,0.4).Figure 12
Robustness simulations of 3-HPA concentration under substrate-sufficient conditions foru=(330,0.4).Figure 13
Robustness simulations of intracellular 1,3-PD concentration under substrate-sufficient conditions foru=(330,0.4).


## 5. Conclusions

Glycerol bioconversion to 1,3-PD byK. pneumoniae in continuous cultures under anaerobic conditions was investigated. Contrasting with the existing models, the paper proposed a new mathematical model by considering 3-HPA inhibitory effects on cells growth and on the activities of the enzymes GDHt and PDOR. Then, we discussed some properties of the system. Furthermore, we presented a parameter identification model to determine the kinetic parameters in the presented system. Since the identification model is subject to constraint of continuous state, we developed a computational approach to solve the parameter identification model based on the penalty function technique and an extension of the state space method. More importantly, a quantitative robustness analysis method was given to infer whether the model is plausible. Numerical results showed that the proposed system is of good robustness.



---

*Source: 690587-2012-11-14.xml*