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

1School of Mathematics and Computer Science, Fujian Normal University, Fuzhou, Fujian 350108, China 2School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan, Hubei 430074, China 3Department of Mathematics and Statistics, Curtin University, Perth, WA 6845, Australia 4School of Mathematical Sciences, Xiamen University, Xiamen, Fujian 361005, China 5College of Life Science, Fujian Normal University, Fuzhou, Fujian 350108, China 6School of Mathematical Sciences, Dalian University of Technology, Dalian, Liaoning 116024, China


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 as Klebsiella pneumoniae, Clostridium butyricum, and Citrobacter freundii have been used to synthesize 1,3-PD from glycerol [1][2][3], among which Klebsiella 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][6][7][8][9][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][13][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][16][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][20][21][22][23][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 engineered K. 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 Section 2, we briefly introduce the metabolic process of glycerol in K. 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.

Glycerol Metabolism in K. pneumoniae and Its Explicit Mathematical Model
During glycerol metabolism by K. 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.
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.,  5 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.
The specific enzyme activities  GDHt ,  PDOR , and  GDH , are expressed as where  is the specific cellular growth rate, and the values of the parameters are listed in Table 2 according to [11,29].
In [11], Sun et al. used the average concentration of the enzymes GDHt, PDOR, and GDH ( protein ) 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]: Here,   max is the maximum concentration of the enzyme ,  = GDHt, PDOR, GDH and  ,  is the inhibitor constant for 3-HPA to the concentration of the corresponding enzyme.

Structural Kinetic Modelling
Let  denote the parameters involved in (3)- (7).The variables   1 and   3 are also regarded as parameters and absorbed into  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  and .In what follows, we will rewrite the vector function  by (, ).The stability of a steady state  0 is determined by the Jacobian matrix of the right hand of (2) with  =  0 , denoted by ( 0 , ).More precisely, the steady state  0 is stable if the largest real part of the eigenvalues of ( 0 , ) is negative, whereas an eigenvalue with a positive real part implies the instability of  0 .
Traditionally, the precise values of the kinetic parameters are required for evaluating the Jacobian ( 0 , ).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 components   of the vector  of orders ≤  + 2 ( is an integer no less than 2) exist and are continuous in  and .
Similarly to [19], taking the following transformation: we can rewrite the system (2) in terms of variables x := ( 1 ,  2 ,  3 )  and rates u(x) = ( 1 (x),  2 (x),  3 (x),  4 (x),  5 (x))  as follows: where Λ = (Λ  ) is a 3 × 5 matrix.The Jacobian of the system (11) at the normalized steady state x = x 0 := (1, 1, 1)  (corresponding to  =  0 ) is where the elements of the matrix Θ u x 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 Θ u x can be expressed as Remark 1.The elements of the matrix Θ u x 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: There are only two independent steady-state reaction rates from the analysis of the stoichiometric matrix .Here, we specify  1 and  2 as the two independent steady-state reaction rates.Let  ∈ (0, 1) denote the percentage of glycerol consumed in the reductive pathway; that is,  :=  2 / 1 .It then follows from the relationship that Finally, setting  1 = V, we can obtain the construction of the matrix Λ as follows: So far, we have provided the parametric presentation of the Jacobian  x .

Specification of the Admissible Space of the Associated Parameters
To evaluate the possible dynamics of the metabolic system on the basis of the Jacobian  x 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  0 , it is obvious that each component of  0 should be no smaller than zero.In addition, (17) requires that each component of  0 cannot be equal to zero.This limitation condition is also reasonable from a biological viewpoint since  0 is a steady state.Conservatively, we set the lower extreme for each component of  0 to be a rather low bound, say, 0.001.As a result,  0 is restricted in Besides, the admissible range of  is defined as  := [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  x ; 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,   3.
When a substance is not only a substrate or product for a reaction, but also has a positive or inhibitory regulation Table 3: The ranges of parameters that have common ranges in the four cases.
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   2  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  2 with respect to the state variable  2 and that of the denominator (note that  2 is a rational function of the state vector ).For the strain without genetic manipulation on the gene of GDHt (including Cases 1 and 4), the reaction rate is given by with So the degree of the numerator of  2 with respect to  2 , denoted by   2 , equals 0 and that of the denominator, denoted by   2 , equals to 2. According to the rule proposed by Steuer et al. [19], we have 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   2  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  2 without product inhibition.So the range of Similarly, we can specify the ranges of   3  2 , which, together with that of   2  2 , is given in Table 4.The range of   5  2 is the same for all the cases, which is given in Table 3.

Dynamics Analysis and Prediction in the Parameter Space
To assess the possible dynamic properties of the system in Case 1, the parameters in (17) and the saturation parameters in ( 13) are uniformly sampled from the corresponding intervals for 1 × 10 6 times.For each tuple ( 0 1, ,  0 2, ,  0 3, ,   , 2 , ) (called random realization),  = 1, 2, . . ., 100000, the Jacobian  x is evaluated and the largest real part of its eigenvalues,  max  x , 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 × 10 6 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.
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   2  2 and   3  2 in maintaining the stability of the metabolic system.For each concerning saturation parameter (  2  2 or   3  2 ),  points ( 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 can be calculated, where  = 5000 is the number of the total samples at each point and    is the number of stable samples at point .Figures 2 and 3 depict the relative fraction of stable models as functions of   2  2 and   3  2 , respectively.The figures illustrate that the probability of obtaining a stable steady state decreases as   2  2 increased (i.e., as the strength of 3-HPA inhibition to  2 weaken) but increases as   3   2 increased.The averaged relative fraction of stable models for the selected  points of   2  2 and   3  2 , which are defined as can be computed, being of 55.877% and 55.713%, respectively.Additionally, the variance of the relative fraction of stable models for the selected  points of   2  2 and   3  2 , which are defined as 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   3  2 is greater than that of   2  2 , which implies that the stability of the system is more sensitive to   3  2 .This conclusion is consistent with the experimental finding that PDOR is more sensitive to the higher level of 3-HPA concentration [32].
Finally, we investigate the influence of the ratio of   the value of  from 0.5 to 2 with a step size of 0.02 and equidistantly select 200 values of   3  2 from [−2, −1] for each point of  (  2  2 will then be uniquely determined by  and   3  2 ).For each pair of (,   3  2 ), randomly generate 5000 values for all other parameters from the corresponding intervals.So we now have 10 6 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  2 to stronger inhibition upon  3 .

Conclusions and Discussions
In this paper, structural kinetic modelling is applied to investigate the stability of glycerol metabolic in the possible engineered K. 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   2  2 and   3  2 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.

Figure 2 :Figure 3 :
Figure 2: The relative fraction of stable models as a function of   2  2 .

Figure 4 :Figure 5 :
Figure 4: The largest real part of the eigenvalues of  x at   3  2 = −0.5 but varied   2  2 .

Figure 8 :
Figure 8: The relative fraction of stable models as a function of .
only the substrate of the reaction   and has no inhibition onto this reaction, and      ∈ [−1, 0] for the product   of the reaction   in the absence of inhibitory feedback from   .Based on this rule, the ranges of the saturation parameters   1  1 ,   2  1 ,   5  1 ,   3  3 ,   4  3 can be given, which are shown in Table

Table 4 :
The ranges of parameters that have different ranges in the four cases.

Table 5 :
Statistical results of the stability of the structural kinetic model.
we calculate the largest real parts of the eigenvalues of the Jacobian  x by varying   2  2 with   3  2 = −0.5 and   3  2 = −1.5, respectively.Similar work is performed on   3  2 with   2  2 fixed.These results are shown in Figures 4-7.Both Figures 4 and 5 illustrate that  max  x is shifted toward positive values as   2 2  2 to   3  2 on the stability of the system.Let  =   2  2 /  3  2 , where   2  2 and   3 2 are both restricted in [−2, −1].We can then obtain that  ∈ [0.5, 2].Then we gradually increase