Abstract

Glycerol can be biologically converted to 1,3-propanediol (1,3-PD) by Klebsiella 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.

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 as Klebsiella pneumoniae, Clostridium butyricum, and Citrobacter freundii have been used to synthesize 1,3-PD from glycerol [13], 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 [510] 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 [1214]. 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 [1517]. 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 [1924]. 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.

2. 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., 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.

Let 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: Then the temporal behavior of the three substances can be described by a set of differential equations: where is a five-dimensional vector function describing the reaction rates of the reactions , and in Figure 1. The components , are defined as follows [11, 29]: The mechanism of (3)–(6) and (7) can be referred to in [11] and [29], respectively. The notations in this work corresponds to in [11]. The physical meanings of the parameters in (3)–(7) are listed as follows:: specific intracellular volume of biomass (L·g−1),: maximum specific transport rate of glycerol (mmol·g−1·h−1),: Michaelis-Menten constant of permease (mmol·L−1),: extracellular concentration of glycerol (mmol·L−1),: with the surface area per biomass (mm2·g−1), the diffusion coefficient of glycerol (L·mm−1·h−1), and the thickness of cell membrane (mm),: conversion coefficient for enzyme activity from μmol·mL−1·min−1 to mmol·L−1·h−1,: average concentration of the enzymes GDHt, PDOR, and GDH in vitro (mg·mL−1),: the ratio of in vivo and in vitro activities of GDHt, PDOR, and GDH,: specific enzyme activities of GDHt, PDOR, and GDH in vitro (U·mg−1),: Michaelis-Menten constants of the enzymes GDHt, PDOR, and GDH (mmol·L−1),: inhibitor constants for 3-HPA to the activities of enzymes GDHt, PDOR and GDH (mmol·L−1),: with the diffusion coefficient of 1,3-PD (L·mm−1·h−1),: 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.

The specific enzyme activities , and , 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 () 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, is the maximum concentration of the enzyme , and is the inhibitor constant for 3-HPA to the concentration of the corresponding enzyme.

3. Structural Kinetic Modelling

Let denote the parameters involved in (3)–(7). The variables and 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 is determined by the Jacobian matrix of the right hand of (2) with , denoted by . More precisely, the steady state is stable if the largest real part of the eigenvalues of is negative, whereas an eigenvalue with a positive real part implies the instability of .

Traditionally, the precise values of the kinetic parameters are required for evaluating the Jacobian . 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 ( 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 and rates as follows: where is a matrix. The Jacobian of the system (11) at the normalized steady state (corresponding to ) is where the elements of the matrix 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 can be expressed as

Remark 1. The elements of the matrix 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 and as the two independent steady-state reaction rates. Let denote the percentage of glycerol consumed in the reductive pathway; that is, . It then follows from the relationship that Finally, setting , we can obtain the construction of the matrix as follows: So far, we have provided the parametric presentation of the Jacobian .

4. Specification of the Admissible Space of the Associated Parameters

To evaluate the possible dynamics of the metabolic system on the basis of the Jacobian 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 , it is obvious that each component of should be no smaller than zero. In addition, (17) requires that each component of cannot be equal to zero. This limitation condition is also reasonable from a biological viewpoint since is a steady state. Conservatively, we set the lower extreme for each component of to be a rather low bound, say, 0.001. As a result, is restricted in Besides, the admissible range of is defined as based on previous research [31]. Furthermore, it is clear that the value of (a positive number) would not influence the sign of the eigenvalues of ; 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, if is only the substrate of the reaction and has no inhibition onto this reaction, and for the product of the reaction in the absence of inhibitory feedback from . Based on this rule, the ranges of the saturation parameters , , , , can be given, which are shown in Table 3.

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 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 with respect to the state variable and that of the denominator (note that 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 with respect to , denoted by , equals and that of the denominator, denoted by , equals to . 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 . Thus, the range of is . 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 without product inhibition. So the range of should be in this context.

Similarly, we can specify the ranges of , which, together with that of , is given in Table 4. The range of is the same for all the cases, which is given in Table 3.

5. 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 times. For each tuple (called random realization), , the Jacobian is evaluated and the largest real part of its eigenvalues, , 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 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 and in maintaining the stability of the metabolic system. For each concerning saturation parameter ( or ), 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 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 and , respectively. The figures illustrate that the probability of obtaining a stable steady state decreases as increased (i.e., as the strength of 3-HPA inhibition to weaken) but increases as increased. The averaged relative fraction of stable models for the selected points of and , 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 and , 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 is greater than that of , which implies that the stability of the system is more sensitive to . This conclusion is consistent with the experimental finding that PDOR is more sensitive to the higher level of 3-HPA concentration [32].

Secondly, given , , , , , , , and , we calculate the largest real parts of the eigenvalues of the Jacobian by varying with and , respectively. Similar work is performed on with fixed. These results are shown in Figures 47. Both Figures 4 and 5 illustrate that is shifted toward positive values as increased whereas Figures 6 and 7 obtain the opposite results for . Furthermore, comparing Figure 4 with Figure 5, we find that the system is stable within a broader interval of while , that is, under weak inhibition to . The opposite conclusion is obtained for 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 at and are, respectively, and , 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 in for the appearance of Hopf bifurcation.

Finally, we investigate the influence of the ratio of to on the stability of the system. Let , where and are both restricted in . We can then obtain that . 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 from for each point of ( will then be uniquely determined by and ). For each pair of , randomly generate 5000 values for all other parameters from the corresponding intervals. So we now have 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 , that is, the system dramatically varies in terms of stability under the transition from stronger inhibition upon to stronger inhibition upon .

6. 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 and 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported in part by National Natural Science Foundation of China (Grant nos. 11301081, 11171050, 11171079, and 31100670), China Postdoctoral Science Foundation (Grant no. 2014M552027), Fujian Department of Education Technology Project (Grant no. JB12041), and Provincial National Science Foundation of Fujian (Grant no. 2014J05001).