# Dynamics Analysis and Prediction of Genetic Regulation in Glycerol Metabolic Network via Structural Kinetic Modelling

Authors: Jianxiong Ye; Honglei Xu; Xueying Huang; Chongrong Ke; Enmin Feng
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Glycerol can be biologically converted to 1,3-propanediol (1,3-PD) byKlebsiella pneumoniae. In the synthesis pathway of 1,3-PD, the accumulation of an intermediary metabolite 3-hydroxypropionaldehyde (3-HPA) would cause an irreversible cessation of the dynamic system. Genetic manipulation on the key enzymes which control the formation rate and consumption rate of 3-HPA would decrease the accumulation of 3-HPA, resulting in nonlinear regulation on the dynamic system. The interest of this work is to focus on analyzing the influence of 3-HPA inhibition on the stability of the dynamic system. Due to the lack of intracellular knowledge, structural kinetic modelling is applied. On the basis of statistical account of the dynamical capabilities of the system in the parameter space, we conclude that, under weak or no inhibition to the reaction of 3-HPA consumption, the system is much easier to obtain a stable state, whereas strong inhibition to its formation is in favor of stabilizing the system. In addition, the existence of Hopf bifurcation in this system is also verified. The obtained results are helpful for deeply understanding the metabolic and genetic regulations of glycerol fermentation by Klebsiella pneumoniae.


---

## Body

## 1. Introduction

1,3-Propanediol (1,3-PD) has wide applications for a variety of markets, especially as a monomer for polyesters, polyethers, and polyurethanes [1]. Microbial production of 1,3-PD, a socially beneficial route to obtain chemicals from renewable resources, has been widely investigated and considered as a competitor to the traditional petrochemical routes. Over the past years, some microorganisms such asKlebsiella pneumoniae,Clostridium butyricum, andCitrobacter freundii have been used to synthesize 1,3-PD from glycerol [1–3], among whichKlebsiella pneumoniae (K. pneumoniae) was most popularly investigated due to its high productivity [4].Over the past years, modelling, stability, and optimal control of glycerol fermentation have been extensively investigated. For details, see [5–10] and the references therein. However, most of the previous researches were based on explicit models for the extracellular substance concentrations instead of a consideration of the metabolic process in the intracellular environment which, in essence, is the origin of multistationarity and oscillation. Up to 2008, the intracellular dynamics of this bioprocess was firstly studied by Sun et al. [11]. Thereafter, we introduced several quantitative measures of biological robustness to estimate the kinetic parameters of Sun’s model in the context where the intracellular data were limited and some of the metabolic mechanisms were not completely known [12–14]. The problem is that huge computing workload will be faced in the evaluation of the biological robustness index, even though the scale of the metabolic network is not too large.In the reductive pathway of glycerol conversion to 1,3-PD, the accumulation of an intermediary metabolite, 3-hydroxypropion-aldehyde (3-HPA), would greatly affect the productivity of the fermentation. There have been many researches on the genetic engineering technology about decreasing 3-HPA accumulation by overexpressing the related genes of the strains [15–17]. However, some questions would emerge on the engineered strains: what kind of changes may happen to the new strains or whether the new strains can inherit “good” properties from the host strains after genetic manipulation. Therefore, it is necessary to reevaluate the dynamics of the engineered strains. In particular, due to the presence of the negative feedback in the 3-HPA accumulation, it is worthwhile to analyze the role of 3-HPA inhibition in the metabolic system of glycerol fermentation.Over the past decades, detailed kinetic models have been widely applied to investigate the dynamic properties of metabolic system, in which differential equations are employed to describe the temporal behavior of the system. Unfortunately, there is a disproportion between the high number of parameters contained in the kinetic models and the relatively incomplete data available [18]. To address these problems, researchers have devoted great effort to exploit other methods to quantitatively evaluate dynamic properties of metabolic processes in the recent years [19–24]. It is worth mentioning that Steuer et al. [19] proposed a structural kinetic modelling approach, which has an advantage in drawing quantitative conclusions about the possible dynamics of the system without assuming detailed knowledge of the underlying enzyme-kinetic rate equations and parameters. One other advantage of this approach is that the classical Michaelis-Menten kinetic model can be transformed to this model scheme with saturation parameters well defined [19].In this paper, our interest is to focus on investigating how 3-HPA inhibition would influence the stability of the glycerol metabolic system of the engineeredK. pneumoniae, in which the genes encoding two key enzymes glycerol dehydratase (GDHt) and 1,3-PD oxydoreductase (PDOR) are overexpressed. A structural kinetic model (SKM) is developed, which is presented in a parametric form with the ranges of all associated parameters well defined. A statistical exploration on the proposed model is performed on the comprehensive parameter space to investigate the system’s capability to obtain a stable state. Additionally, by varying the strength of 3-HPA inhibitions upon its upstream and downstream reactions, we verify the existence of Hopf bifurcation. Finally, we analyze how the ratio of GDHt activity to PDOR activity may influence the stability of the system.This paper is organized as follows. In Section2, we briefly introduce the metabolic process of glycerol inK. pneumoniae and present the explicit rate equations for the considered metabolites. In Section 3, a SKM is developed for the reductive pathway of this process. In Section 4, the physiologically feasible ranges of the associated parameters in the SKM are specified. In Section 5, the stability of the system is statistically investigated on the comprehensive parameter space. Conclusions and discussions of the computational results are presented at the end of this paper.


## 2. Glycerol Metabolism inK. pneumoniae and Its Explicit Mathematical Model

During glycerol metabolism byK. pneumoniae under anaerobic condition, glycerol is first transported across cell membrane from the extracellular environment to the intracellular environment. In the intracellular environment, glycerol is dissimilated through coupled oxidative and reductive pathways as shown in Figure 1. The goal product 1,3-PD is produced by the reductive branch in two successive enzymatic reactions [25]: glycerol is first dehydrated to 3-HPA by the enzyme GDHt; 3-HPA is then converted to 1,3-PD by the enzyme PDOR. In the oxidative pathway, glycerol oxidation is catalyzed by the enzyme glycerol dehydrogenase (GDH), leading to the formation of dihydroxyacetone (DHA); DHA is further phosphorylated by two dihydroxyacetone kinases and is channelled into glycolysis, yielding the same fermentation products as in glucose fermentation (acetate, ethanol) and to the generation of energy and reducing power. Finally, the products are transported across cell membrane from the intracellular environment to the extracellular environment.Figure 1
Pathway of anaerobic glycerol metabolism inK. pneumoniae. Solid line: metabolic reaction; dashed line: catalytic action or inhibitory regulation; dashed-dot line: the omitted intermediary metabolic reactions; ⊕: positive regulation; ⊖: inhibitory regulation.The flux through the reductive pathway (i.e., 1,3-PD synthesis pathway) is largely controlled by the synthesis of the enzymes GDHt and PDOR, because the accumulation of 3-HPA represses the expression of the genes of the enzymes GDHt and PDOR as well as that of the enzyme GDH [26]. In addition, the accumulation of 3-HPA also has inhibitory effects on the activities of the above three enzymes [25, 27, 28]. Therefore, there appears a negative feedback in 3-HPA accumulation.The accumulation of 3-HPA in the reductive pathway can be decreased by overexpressing the genes of the two enzymes in this branch (i.e., the enzymes GDHt and PDOR) [15]. According to the genetic manipulation on the two enzymes, the strains discussed in this paper can be divided into four cases.Case 1.
The original strain (no genetic manipulation).Case 2.
In the constructed strain, both GDHt and PDOR are overexpressed.Case 3.
In the constructed strain, only GDHt is overexpressed.Case 4.
In the constructed strain, only PDOR is overexpressed.As we mentioned above, the goal product 1,3-PD is synthesized in the reductive pathway and the accumulation of 3-HPA in this branch would greatly affect the productivity of the fermentation and the stability of the system. Our main concern in this work is therefore focused on the reductive pathway. In addition, the oxidative reaction of glycerol (i.e.,r5 in Figure 1) is also considered. One reason for introducing this reaction is to guarantee the material balance. Another reason is that the accumulation of 3-HPA also inhibits this reaction, which would affect the stability of the system.LetS∶=(S1,S2,S3)T denote the intracellular concentrations of glycerol, 3-HPA and 1,3-PD. According to Figure 1, we can obtain the stoichiometric matrix covering the formations and consumptions of glycerol, 3-HPA and 1,3-PD as follows:
(1)N=1-100-101-100001-10.
Then the temporal behavior of the three substances can be described by a set of differential equations:
(2)dSdt=N·V,
where V=(V1,V2,V3,V4,V5)T is a five-dimensional vector function describing the reaction rates of the reactions r1,r2,r3,r4, and r5 in Figure 1. The components Vj, j=1,2,3,4,5 are defined as follows [11, 29]:
(3)V1=1VsJmax⁡S1eS1e+Km+1AsS1e-S1,(4)V2=bCproteinKGDHtUGDHtS1KmGDHt(1+S2/KiGDHt)+S1,(5)V3=bCproteinKPDORUPDORS2KmPDOR+S2(1+S2/KiPDOR),(6)V4=KPDS3-S3e,(7)V5=bCproteinUGDHS1(1+S2/KiGDH)(KmGDH+S1).
The mechanism of (3)–(6) and (7) can be referred to in [11] and [29], respectively. The notations S1,S2,S3 in this work corresponds to Csi,C3-HPA,CPDi in [11]. The physical meanings of the parameters in (3)–(7) are listed as follows:Vs: specific intracellular volume of biomass (L·g−1),Jmax⁡: maximum specific transport rate of glycerol (mmol·g−1·h−1),Km: Michaelis-Menten constant of permease (mmol·L−1),S1e: extracellular concentration of glycerol (mmol·L−1),As: As∶=δ/B·Df with B the surface area per biomass (mm2·g−1), Df the diffusion coefficient of glycerol (L·mm−1·h−1), and δ the thickness of cell membrane (mm),b: conversion coefficient for enzyme activity from μmol·mL−1·min−1 to mmol·L−1·h−1,Cprotein: average concentration of the enzymes GDHt, PDOR, and GDHin vitro (mg·mL−1),KGDHt,KPDOR,KGDH: the ratio ofin vivo andin vitro activities of GDHt, PDOR, and GDH,UGDHt,UPDOR,UGDH: specific enzyme activities of GDHt, PDOR, and GDHin vitro (U·mg−1),KmGDHt,KmPDOR,KmGDH: Michaelis-Menten constants of the enzymes GDHt, PDOR, and GDH (mmol·L−1),KiGDHt,KiPDOR,KiGDH: inhibitor constants for 3-HPA to the activities of enzymes GDHt, PDOR and GDH (mmol·L−1),KPD: KPD∶=B·Df′/Vs·δ with Df′ the diffusion coefficient of 1,3-PD (L·mm−1·h−1),S3e: extracellular concentration of 1,3-PD (mmol·L−1).The values of the above parameters are roughly estimated in [11, 29], which are listed in Table 1.Table 1
The values of parameters in (3)–(7).
V
s
J
max
⁡
K
m
A
s
b
C
protein
0.151
54.664
1.34
1.896
×
1
0
-
4
60
1.7
K
GDHt
K
PDOR
K
GDH
K
m
GDHt
K
m
PDOR
K
m
GDH
0.171
0.217
0.274
0.53
0.14
1.27
K
i
GDHt
K
i
PDOR
K
i
GDH
K
PD
—
—
220.319
0.418
986.251
25.137
—
—The specific enzyme activitiesUGDHt,UPDOR, and UGDH, are expressed as
(8)Ui=Ui0+aiμ+ΔUimS1eS1e+KiS,i=GDHt,PDOR,GDH,
where μ is the specific cellular growth rate, and the values of the parameters are listed in Table 2 according to [11, 29].Table 2
The values of parameters in (8).
i
U
i
0
a
i
Δ
U
i
m
K
i
S
GDHt
2.8
-
3.53
-
0.87
2.41
PDOR
2.32
−1.87
−1.38
1.35
GDH
14.05
−3.94
−9.76
2.40In [11], Sun et al. used the average concentration of the enzymes GDHt, PDOR, and GDH (Cprotein) to approximate their actual concentrations. In consideration of the inhibitory effects of 3-HPA accumulation onto the DNA synthesis of the three enzymes [15], we take the concentration of each enzyme as a function of 3-HPA, which is given as follows [26]:
(9)Cproteini=Cmax⁡i1+S2/Kmp,i,i=GDHt,PDOR,GDH.
Here, Cmax⁡i is the maximum concentration of the enzyme i, i=GDHt,PDOR,GDH and Kmp,i is the inhibitor constant for 3-HPA to the concentration of the corresponding enzyme.


## 3. Structural Kinetic Modelling

LetK denote the parameters involved in (3)–(7). The variables S1e and S3e are also regarded as parameters and absorbed into K due to the characteristic features of the structural kinetic modelling, which will be shown in the next section.Observing the right-hand side of (3)–(7), we can find that the reaction rates are represented by rational functions of S and K. In what follows, we will rewrite the vector function V by V(S,K). The stability of a steady state S0 is determined by the Jacobian matrix of the right hand of (2) with S=S0, denoted by J(S0,K). More precisely, the steady state S0 is stable if the largest real part of the eigenvalues of J(S0,K) is negative, whereas an eigenvalue with a positive real part implies the instability of S0.Traditionally, the precise values of the kinetic parameters are required for evaluating the JacobianJ(S0,K). Although the parameters involved in (3)–(7) have been estimated based on 58 groups of steady-state experimental data [11], it is hard to verify the reliability of these parameter values due to the lack of process data of the intracellular substances. In fact, the process data of the intracellular substances cannot be measured using the existing technology. Furthermore, the parameters values may depend on many factors such as strain types or experimental and physiological conditions. In other words, the obtained parameter values from the literature may be invalid for the new strains. To overcome this problem, we will use the method of structural kinetic modelling proposed by Steuer et al. [19] to investigate the stability of the metabolic system, because this method can be used to draw qualitative conclusions about the possible dynamics of the system without very detailed knowledge of the underlying enzyme-kinetic rate equations and parameters.We make the following assumption.(H1)
All partial derivatives of the componentsVj of the vector V of orders ≤L+2 (L is an integer no less than 2) exist and are continuous in S and K.Similarly to [19], taking the following transformation:
(10)xi∶=SiSi0,Λij=NijVj(S0,K)Si0,μjx∶=Vj(S,K)Vj(S0,K),i=1,2,3,j=1,2,3,4,5
we can rewrite the system (2) in terms of variables x∶=(x1,x2,x3)T and rates u(x)=(μ1(x),μ2(x),μ3(x),μ4(x),μ5(x))T as follows:
(11)dxdt=Λux,
where Λ=(Λij) is a 3×5 matrix. The Jacobian of the system (11) at the normalized steady state x=x0∶=(1,1,1)T (corresponding to S=S0) is
(12)Jx=Λ∂u(x)∂xx=x0∶=ΛΘxu,
where the elements of the matrix Θxu specify the normalized saturation degrees of the reactions with respect to their substrates or products, defined as saturation parameters [20]. According to the metabolic network in Figure 1, the parametric presentation of the matrix Θxu can be expressed as
(13)Θxu=θx1μ100θx1μ2θx2μ200θx2μ3θx3μ300θx3μ4θx1μ5θx2μ50.Remark 1.
The elements of the matrixΘxu specify the normalized degree of saturation of each reaction with respect to each of the reactants, defined in closed analogy to the scaled elasticity coefficients of metabolic control analysis. In other words, each of them describes the degree of the influence of the substrate concentration on the specified reaction rate. Therefore, if the substrate has strong regulation onto the specified reaction, the absolute value of the corresponding saturation parameter should be large; if the substrate has weak regulation onto the specified reaction, the absolute value of the corresponding parameter should be small.At steady state, none of the state variables changes in value, even though material is flowing through the system. Mathematically, this situation requires that all derivatives are zero:(14)dSdt=N·V=0.
There are only two independent steady-state reaction rates from the analysis of the stoichiometric matrix N. Here, we specify V1 and V2 as the two independent steady-state reaction rates. Let α∈(0,1) denote the percentage of glycerol consumed in the reductive pathway; that is, α∶=V2/V1. It then follows from the relationship
(15)V2+V5=V1
that
(16)V5=1-αV1.
Finally, setting V1=v, we can obtain the construction of the matrix Λ as follows:
(17)Λ=vS10-αvS1000-(1-α)vS100αvS20-αvS200000αvS30-αvS300.
So far, we have provided the parametric presentation of the Jacobian Jx.


## 4. Specification of the Admissible Space of the Associated Parameters

To evaluate the possible dynamics of the metabolic system on the basis of the JacobianJx in (12), we should firstly specify the feasible ranges of the associated parameters. According to the practical experiments, the critical concentration of glycerol, 3-HPA and 1,3-PD are 2039, 30, and 939.5 mmol/L [2, 30], respectively. As for the lower bound of S0, it is obvious that each component of S0 should be no smaller than zero. In addition, (17) requires that each component of S0 cannot be equal to zero. This limitation condition is also reasonable from a biological viewpoint since S0 is a steady state. Conservatively, we set the lower extreme for each component of S0 to be a rather low bound, say, 0.001. As a result, S0 is restricted in
(18)Wad∶=0.001,2039×0.001,30×0.001,939.5.
Besides, the admissible range of α is defined as I∶=[0.55,0.65] based on previous research [31]. Furthermore, it is clear that the value of v (a positive number) would not influence the sign of the eigenvalues of Jx; that is, it is irrelative to the stability of the system. So we assign its value to be 100 mmol·L−1·h−1 in the remaining of this paper.The ranges of saturation parameters in (13) can be specified by the rule proposed by Steuer et al. [19]. In particular, θxiμj∈[0,1] if xi is only the substrate of the reaction rj and has no inhibition onto this reaction, and θxiμj∈[-1,0] for the product xi of the reaction rj in the absence of inhibitory feedback from xi. Based on this rule, the ranges of the saturation parameters θx1μ1, θx1μ2, θx1μ5, θx3μ3, θx3μ4 can be given, which are shown in Table 3.Table 3
The ranges of parameters that have common ranges in the four cases.
θ
x
1
μ
1
θ
x
1
μ
2
θ
x
3
μ
3
θ
x
3
μ
4
θ
x
1
μ
5
θ
x
2
μ
5
[
-
1,0
]
[
0,1
]
[
-
1,0
]
[
0,1
]
[
0,1
]
[
-
2,0
]When a substance is not only a substrate or product for a reaction, but also has a positive or inhibitory regulation on that reaction, the range of the corresponding saturation parameter needs to be determined by the explicit expression of the rate equation for that reaction. We take the parameterθx2μ2 as an example to illustrate how to define the range of the saturation parameter of this type. The required knowledge is only the degree of the numerator of V2 with respect to the state variable S2 and that of the denominator (note that V2 is a rational function of the state vector S). For the strain without genetic manipulation on the gene of GDHt (including Cases 1 and 4), the reaction rate is given by
(19)V2=bCproteinGDHtKGDHtUGDHtS1KmGDHt(1+S2/KiGDHt)+S1
with
(20)CproteinGDHt=Cmax⁡GDHt1+S2/Kmp,GDHt.
So the degree of the numerator of V2 with respect to S2, denoted by nV2, equals 0 and that of the denominator, denoted by mV2, equals to 2. According to the rule proposed by Steuer et al. [19], we have
(21)θx2μ2∶=nV2-γmV2,γ∈0,1,
where γ corresponds to the notation α that appeared in (8) of the supporting appendix of [19]. It was proved in this supporting appendix that γ∈[0,1]. Thus, the range of θx2μ2 is [-2,0]. Alternatively, for the strain with genetic manipulation on the gene of GDHt (including Cases 2 and 3), the inhibitory effect of 3-HPA accumulation onto the synthesis of GDHt is eliminated. In other words, 3-HPA is only the product of the reaction r2 without product inhibition. So the range of θx2μ2 should be [-1,0] in this context.Similarly, we can specify the ranges ofθx2μ3, which, together with that of θx2μ2, is given in Table 4. The range of θx2μ5 is the same for all the cases, which is given in Table 3.Table 4
The ranges of parameters that have different ranges in the four cases.
Case1
Case2
Case3
Case4
θ
x
2
μ
2
[
-
2,0
]
[
-
1,0
]
[
-
1,0
]
[
-
2,0
]
θ
x
2
μ
3
[
-
2,1
]
[
0,1
]
[
-
2,1
]
[
0,1
]


## 5. Dynamics Analysis and Prediction in the Parameter Space

To assess the possible dynamic properties of the system in Case1, the parameters in (17) and the saturation parameters in (13) are uniformly sampled from the corresponding intervals for 1×106 times. For each tuple (S1,k0,S2,k0,S3,k0,αk,θx1,kμ1,θx1,kμ2,θx2,kμ2,θx2,kμ3,θx3,kμ3,θx3,kμ4,θx1,kμ5,θx2,kμ5) (called random realization), k=1,2,…,100000, the Jacobian Jx is evaluated and the largest real part of its eigenvalues, λJxmax⁡, is recorded. Consequently, the stability of the system can be obtained for each random realization based on the classical theory of dynamical systems. The percentage of stable models among the 1×106 random realizations is then accounted. Likewise, we statistically evaluate the percentage of stable models for the other three cases and the results are listed in Table 5.Table 5
Statistical results of the stability of the structural kinetic model.
Stability
Case1
Case2
Case3
Case4
Stable
56.2%
85.2%
38.8%
91.5%
Unstable
43.8%
14.8%
61.2%
8.5%To investigate the role of 3-HPA inhibition in more detail and verify the existence of Hopf bifurcation, we carry out the following analysis.Firstly, we, respectively, discuss the roles ofθx2μ2 and θx2μ3 in maintaining the stability of the metabolic system. For each concerning saturation parameter (θx2μ2 or θx2μ3), Q points (Q is set to be 200 in the actual numerical implementation) are equidistantly selected from its interval with the other parameters sampled from a uniform distribution (5000 samples in total for each selected point of the concerning parameter). Then the relative fraction of stable models, which is defined as
(22)pij=MijM,i=θx2μ2,θx2μ3,j=1,2,…,Q,
can be calculated, where M=5000 is the number of the total samples at each point and Mij is the number of stable samples at point j. Figures 2 and 3 depict the relative fraction of stable models as functions of θx2μ2 and θx2μ3, respectively. The figures illustrate that the probability of obtaining a stable steady state decreases as θx2μ2 increased (i.e., as the strength of 3-HPA inhibition to r2 weaken) but increases as θx2μ3 increased. The averaged relative fraction of stable models for the selected Q points of θx2μ2 and θx2μ3, which are defined as
(23)u¯i=∑j=1Q‍pijQ,i=θx2μ2,θx2μ3,
can be computed, being of 55.877% and 55.713%, respectively. Additionally, the variance of the relative fraction of stable models for the selected Q points of θx2μ2 and θx2μ3, which are defined as
(24)σi=∑j=1Q‍pij-u¯i2Q-1,i=θx2μ2,θx2μ3,
can also be computed, which are 0.041 and 0.108, respectively. The results indicate that the variance of the relative fraction of stable models caused by the change of θx2μ3 is greater than that of θx2μ2, which implies that the stability of the system is more sensitive to θx2μ3. This conclusion is consistent with the experimental finding that PDOR is more sensitive to the higher level of 3-HPA concentration [32].Figure 2
The relative fraction of stable models as a function ofθx2μ2.Figure 3
The relative fraction of stable models as a function ofθx2μ3.Secondly, givenS0=(400,20,200)T, θx1μ1=-0.8, θx1μ2=0.5, θx3μ3=-0.8, θx3μ4=0.8, θx1μ5=0.5, θx2μ5=-1.0, and α=0.6, we calculate the largest real parts of the eigenvalues of the Jacobian Jx by varying θx2μ2 with θx2μ3=-0.5 and θx2μ3=-1.5, respectively. Similar work is performed on θx2μ3 with θx2μ2 fixed. These results are shown in Figures 4–7. Both Figures 4 and 5 illustrate that λJxmax⁡ is shifted toward positive values as θx2μ2 increased whereas Figures 6 and 7 obtain the opposite results for θx2μ3. Furthermore, comparing Figure 4 with Figure 5, we find that the system is stable within a broader interval of θx2μ2 while θx2μ3=-0.5, that is, under weak inhibition to r3. The opposite conclusion is obtained for θx2μ3 based on the comparison between Figures 6 and 7. Additionally, in the numerical computation of Figure 5, we obtain that the eigenvalues of the Jacobian Jx at θx2μ2=-1.32 and θx2μ2=-1.30 are, respectively, {-0.00927539+0.776809i,-0.00927539-0.776809i,-0.29554} and {0.0183351+0.761386i,0.0183351-0.761386i,-0.296216}, which shows that the sign of the real part of the pair of complex conjugate eigenvalues has been changed. Under assumption (H1) and applying Hopf bifurcation theorem [33], we conclude that there exists a critical value θ¯x2μ2 in (-1.32,-1.30) for the appearance of Hopf bifurcation.Figure 4
The largest real part of the eigenvalues ofJx at θx2μ3=-0.5 but varied θx2μ2.Figure 5
The largest real part of the eigenvalues ofJx at θx2μ3=-1.5 but varied θx2μ2.Figure 6
The largest real part of the eigenvalues ofJx at θx2μ2=-0.5 but varied θx2μ3.Figure 7
The largest real part of the eigenvalues ofJx at θx2μ2=-1.5 but varied θx2μ3.Finally, we investigate the influence of the ratio ofθx2μ2 to θx2μ3 on the stability of the system. Let β=θx2μ2/θx2μ3, where θx2μ2 and θx2μ3 are both restricted in [-2,-1]. We can then obtain that β∈[0.5,2]. Then we gradually increase the value of β from 0.5 to 2 with a step size of 0.02 and equidistantly select 200 values of θx2μ3 from [-2,-1] for each point of β (θx2μ2 will then be uniquely determined by β and θx2μ3). For each pair of (β,θx2μ3), randomly generate 5000 values for all other parameters from the corresponding intervals. So we now have 106 realizations at each point of β. After that, we account the percentage of stable models among the sampled realizations for each β. As shown in Figure 8, the percentage of stable models is monotonously increasing with respect to β. Furthermore, the curve shows that the probability of finding a stable state is markedly increased in the neighborhood of β=1, that is, the system dramatically varies in terms of stability under the transition from stronger inhibition upon r2 to stronger inhibition upon r3.Figure 8
The relative fraction of stable models as a function ofβ.


## 6. Conclusions and Discussions

In this paper, structural kinetic modelling is applied to investigate the stability of glycerol metabolic in the possible engineeredK. pneumoniae. On the basis of the numerical results, we can summarize some essential properties of the reductive pathway. Firstly, both the results in Table 5 and Figure 8 reveal that, under stronger 3-HPA inhibition to its formation than that to its consumption, the system is much easier to obtain a stable state. Secondly, by varying θx2μ2 and θx2μ3 separately, we find that the stability of the system is much more sensitive to the latter, which confirms the previous experimental research that PDOR is more sensitive to higher level of 3-HPA concentration [11, 32]. Thirdly, our numerical result shows the existence of Hopf bifurcation, which is also consistent with the previous theoretical work of the extracellular reactants [7]. Finally, combining the previous experimental researches with our theoretical analysis, we conjecture that overexpression of the genes encoding GDHt and PDOR such that the activity of GDHt is higher than that of PDOR is favored for stabilizing the metabolic system, resulting in increasing the productivity of 1,3-PD. The contributions of this work would be helpful for deeply understanding metabolic and genetic regulation of glycerol metabolism.This is so far the first trail to investigate the stability of glycerol metabolic system in the intracellular environment. Although the present work only considers a relatively simple part of the whole network of glycerol metabolism, the obtained results have shown the validity of structural kinetic modelling in investigating the dynamics of the metabolic system. Most important of all, the method needs far less computation cost compared with the techniques in [12, 13]. Therefore, if the temporal changes of the metabolites are not the primary concern, we believe that structural kinetic modelling would be much easier to be extended to study the dynamics of the whole network of glycerol metabolism, which is also one of the emphases of our future research.



---

*Source: 673120-2015-03-02.xml*