Mathematical study of feedback inhibition effects on the dynamics of metabolites on the central metabolism of a yeast cell: a combination of kinetic model and metabolic control analysis

Abstract The fermentation system is a metabolic system that can be considered as a complex system, since it involves metabolites linked by different reactions. Some outputs of the reaction act as input in the other reactions and they control the behaviour of the system. Glucose as the main carbon source of the yeast cells plays an important role in the cell’s reproduction cycle and product synthesis. Glucose promotes negative feedback to some reactions in the central metabolism of the yeast cell. In this paper we mathematically study the effects of the negative feedback (inhibition) by glucose and other metabolites on the dynamics of the fermentation system. We also study the sensitivity of the flux in the presence of inhibition. We found that high inhibition by glucose affects the concentration of acetyl-CoA, which will lead to the respiration pathway of the yeast cells. High inhibition by acetaldehyde affects the concentration of all metabolites, including the maximal concentration of ethanol. Based on the metabolic control analysis results, we found that glucose can be considered as the external regulating point in increasing the flux of ethanol. Although glucose acts as a negative feedback, it can also be used to promote certain processes in the metabolic pathway since it gives the highest positive control in increasing the flux of ethanol as the desired product. For the internal regulation point, the maximal concentration of ethanol can be increased significantly by regulating the maximal activity of alcohol dehydrogenase and acetaldehyde dehydrogenase simultaneously.


Introduction
Introducing a new cellular process is one of the genetic engineering targets for improving cell productivity [1,2]. Metabolic engineering is a process that can convert, for example, renewable carbon sources into biofuels [3]. The process alters enzyme activity to redirect carbon fluxes toward biofuel production [3,4]. Since conversion of substrate into biofuels involves a complex metabolic network, theoretical approach is needed as a tool to study the complex system to complete the experimental study. Theoretical studies provide high predictive power to design cells with desirable properties cheaply and reliably according to the needs of biotechnology production [5,6]. Some researchers used advanced mathematical model to capture the dynamical behaviour of a complex metabolic system such as the fermentation system [7][8][9][10][11][12][13][14]. The model was applied as a predictive tool that can be compared with experimental data in terms of conversion of substrate becomes product under steady state condition. A valid model would presents a good agreement with the experimental data such that the model can be used as a miniature or a model prediction of the complex metabolic system. The model can be used as a strategic tool for the design and optimisation of metabolic network in metabolic engineering [1,8].
Each complex metabolic system is regulated by several enzymes which control the rate of the reaction process. In the biochemical pathway, the product of one reaction becomes a substrate in the next reaction until the final product is synthesised. Fermentation system is a complex metabolic system where yeast, such as Saccharomyces cerevisiae, converts glucose as the substrate to pyruvate acid via the glycolysis pathway and then go one step further converting pyruvate into ethanol. One glucose molecule is converted into two ethanol molecules and two carbon dioxide molecules [15,16]. Each reaction process is catalysed by different enzyme with different mechanism. In the previous study, we derived a kinetic model of the fermentation system by taking into consideration the central metabolism of fermentation pathway. The characteristics of enzymes in the main yeast cell pathway were studied, such as the possibility of delays in a reaction due to the existence of a rate limiting step, in which the conformational changes of the enzyme structure require certain time to rearrange the structure until it is ready to bind the substrate again [17]. We found that there existed a bifurcation diagram for the delay that depended on the rate of glucose supply and kinetic parameters of the delayed reaction. There also existed a certain glucose supply that may optimise ethanol production. In the other study, we also investigated the key enzymes that can be regulated to enhance ethanol production [9]. We found that there were supply glucose as an external control parameter and three enzymes as internal regulation targets that can be adjusted in enhancing ethanol production. The three enzymes are acetaldehyde dehydrogenase, pyruvate decarboxylase and alcohol dehydrogenase which reside in the acetaldehyde branch point. In this study, we will continue our investigation regarding analysis of the optimality conditions of the yeast cell to produce ethanol with high concentration. The existence of internal feedback inhibitions becomes the key point of this study since some experimental researchers found that acetaldehyde acts as inhibitor for glycolytic pathway and glucose inhibits the reaction catalysed by pyruvate dehydrogenase complex [18][19][20]. Here, we continue our theoretical analysis by using mathematical model to study the effects of the existence of inhibitions to the dynamical behaviour of metabolites in the central metabolism of the yeast cell, especially its effects to the optimal production of ethanol. We also apply metabolic control analysis to identify the sensitivity of the positive steady-state properties of the ethanol concentration and ethanol flux to a small parameter change due to the existence of feedback inhibition. Our analysis serves as a starting point to investigate the major biotechnological questions such as the observation of complex cellular behaviours, the prediction of the outcomes of metabolic with perturbations and the improvement of the desired property of the yeast metabolism.
The paper is organized as follows: in the Methods section, we formulate the kinetic model of fermentation system by taking into consideration the internal feedback inhibition by glucose and acetaldehyde and also present the mathematical analysis of the model. In Results and discussion section, some numerical simulations are presented to study the transient behavior of the metabolites and using metabolic control analysis, the steady-state properties of the system are discussed to analyze the sensitivity of kinetic parameters. Some remarks and conclusions are presented in the last section.

Kinetic formulation
In this paper, we propose a mathematical model describing the chemical conversion of metabolites in the central metabolism of a yeast cell. The modelling formulation follows the line of Lei et al. [7] and Kasbawati et al. [9]. We consider a continuous fermentation process that occurred at ideal conditions. As shown in Figure 1, we assume that the carbon metabolism on the central yeast metabolic pathway is started with the uptake of glucose. The glucose (S 1 ) is converted into pyruvate (S 2 ) via glycolysis pathway, a lumped reaction [7]. Chemical conversion of glucose into pyruvate as the last product of glycolytic pathway is modelled using Michaelis-Menten kinetic equation [7,9]. Since there exists uncompetitive internal inhibition by acetaldehyde (S 4 ) to the glycolytic pathway [20], the Michaelis-Menten kinetic equation becomes where K a is an inhibition constant of S 4 with respect to r 1 ; V 1 is the maximum reaction of glucose conversion to become pyruvate in the glycolysis pathway and K 1 is the Michaelis constant. The pyruvate is then converted to become acetaldehyde ðS 4 ) by pyruvatedecarboxylase (r 2 ) and acetyl-CoA (S 3 Þ by pyruvate dehydrogenase complex (r 3 ). The kinetics for pyruvate decarboxylase is with maximum reaction rate V 2 and Michaelis constants K 2 : Based on the study of Alexander and Jeffries [19], it was known that glucose is an inhibitor for the pyruvate dehydrogenase complex. It inhibits the enzyme noncompetitively meaning that the higher the concentration of the glucose, the more likely that glucose will bind to the allosteric site of the enzyme. When the levels of glucose fall below a threshold, the effect of the negative feedback diminishes and pyruvate dehydrogenase complex is reactivated. Therefore, we have the kinetic equation for this reaction: with maximum reaction rate V 3 ; Michaelis constants K 2 ; and inhibition constant K b : The product of r 3 then leads to the respiration pathway of the yeast cell, i.e. the TCA cycle, which produces energy for the cell, whereas acetaldehyde leads to the fermentative pathway, in which it is converted to become acetate (S 5 ) and ethanol (S 6 ). Both conversion processes are, respectively, catalysed by acetaldehyde dehydrogenase and alcohol dehydrogenase. Alcohol dehydrogenase catalyses the reaction by a reversible mechanism, whereas acetaldehyde dehydrogenase catalyses the reaction by an irreversible mechanism, such that we have the kinetic equation for the both reactions: Acetyl-CoA syntehase then converts acetate to become acetyl-CoA that also leads to the respiration pathway. Since it is known that there also exists noncompetitive internal inhibition by glucose to this reaction [21], we have the kinetic equation for this reaction: where K c is an inhibition constant of S 1 with respect to r 6 ; V 6 is the maximum reaction rate of acetyl-CoA syntehase, and K 6 is the Michaelis constant.
By using the mass action rate law and the formulated kinetic equation (1)-(6), we then derive the mathematical model of the metabolites in the central metabolism of the yeast cell with internal inhibition as follows: with initial conditions where S a ;. . . ; S f are initial concentration of metabolites which are greater or equal to zero. All parameters in model (7) are assumed positive. The mathematical properties of the model are derived in the following section.

Mathematical analysis
In this section we use stability theory to derive the mathematical properties of the kinetic model (7). One of the mathematical properties of the system that will be analysed is the stability conditions of the positive steady state of system (7). The steady state solution is defined in the region such that we get the equilibrium solution for S 1 depending on the equilibrium solution of S 4 , The equilibrium solution of S 1 will be positive if it fulfils the conditions These conditions imply that the maximum rate of reaction in the glycolysis pathway should be greater than the maximum inflow of glucose supply, and steady state concentration of S 4 should be maintained in the interval 0; V 1 ÀaG which can be rewritten in quadratic polynomial, The roots of polynomial (12) are equilibrium solutions of S 2 which depend on the equilibrium solutions of S 1 and S 4 . If the equilibrium solution of S 1 fulfils condition (10), then we have a 3 >0. Using Descartes rule of sign [22], we have two possibilities number of positive roots of (12). When a 1 > 0 and a 2 < 0; two positive roots of polynomial (12) exist. Otherwise, when a 1 < 0 and a 2 > 0 or a 2 < 0; one positive root of polynomial (12) exists. Therefore, polynomial (12) will produce at most two positive roots for S 2 . This condition leads to the possibility of the existence of multiple steady states in the system where the fermentation system will generate different steady state concentrations for different initial conditions. Furthermore, from (7c) we get the equilibrium solution of S 3 : that will be positive if the equilibrium solutions of S 1 , S 2 and S 5 are positive. From (7d), we have that can be rewritten as a quadratic polynomial as follows: If the equilibrium solutions of S 2 and S 6 are positive, then we get b 3 >0. Using the similar analysis with S 2 , we find that the polynomial (15) will produce at most two positive roots for S 4 . This result also leads to the possibility of multiple steady states in the modelled fermentation system. Furthermore, from (7e) we have that also can be rewritten in a quadratic polynomial form, i.e. with Since c 1 >0 and c 3 <0 when S 1 and S 4 are positive, then polynomial (17) will produce at most one positive equilibrium solution for S 5 . Furthermore, from the last Equation (7f), we have In quadratic polynomial form, Equation (18) can be rewritten as follows: with Since d 1 >0 and d 3 <0 when S 4 >0, we get that S 4 has only one positive equilibrium solution. Summarising our analysis results, we found that the fermentation system (7) has an equilibrium solution respectively, stand for the i-th positive root of S 2 and S 4 . This steady state solution shows the steady state concentration of metabolites when the inhibitions exist in the system. Furthermore, the mathematical property of the system (7) that will be studied is the local stability of the equilibrium point E. We analyse its local stability behaviour by investigating the eigenvalues of the Jacobian matrix of the linearised system (7). By linearising system (7) around the equilibrium point E, we get the following Jacobian matrix: where We can observe that two eigenvalues of matrix J are negative, i.e.
and k 2 ¼ Àa: The rest of eigenvalues of matrix J is the solutions of the characteristic equation, where coefficients a 1 ; . . . ; a 4 depend on D 1 ; . . . ; D 12 : Because the coefficients of polynomial (22) depend on the steady state E and the large number of kinetic parameters, it is difficult to determine the value of its roots. Routh-Hurwitz criteria gives conditions on the coefficients of a polynomial such that the zeros of the polynomial have negative real part [22]. For polynomial (22), the conditions are i. a 1 > 0; ii. a 1 1 0 a 2 > 0 or a 1 a 2 > 0; iii. Àa 1 a 4 a 2 2 À a 2 a 4 a 3 À a 1 a 4 a 3 > 0: Therefore, the steady state E will be locally asymptotically stable, if it fulfils the conditions (i)-(iv) meaning that for a certain experimental condition, the fermentation system will produce metabolites with the amount of concentration indicated by point E. In the next section we will present some numerical simulations to clarify these mathematical results and to capture the transient behaviour of the system (7) for a certain experimental condition indicated by kinetic values given in Table 1.

Numerical simulations
In this section, we present some numerical simulations to validate the mathematical analysis and to capture the transient behaviour of solutions of system (7). Some simulations are presented to observe effects of the existence of the internal inhibitions. Using kinetic parameter values given in Table 1, we get one positive steady state,  (23) with eigenvalues k 1 ¼ À0:01; k 2 ¼ À0:01; k 3 ¼ À1:22; k 4 ¼ À1:103Â10 7 ; k 5 ¼ À17142:3; k 6 ¼ À0:22 implying the local stability of E. By choosing initial conditions S 1 (0)¼20 g L À1 and zero for the others and final time observation 100 hours, we get numerical solutions of system (7) depicted in Figure 2. Figure 2 shows the transient behaviour of metabolites when the system experiences inhibitions by glucose and acetaldehyde. We can observe that inhibitions of some reactions highly affect the concentration of acetyl-CoA, which will lead to the respiration pathway of the yeast cell. The feedback inhibitions extremely decrease the concentration of acetyl-CoA. However, the feedback inhibition by acetaldehyde to the glycolysis pathway did not influence the production of pyruvate since there is no gap between the solution with and without inhibition term. This result implies that the inhibition by glucose to the third and the sixth reactions should be taken into attention since it influences the respiratory pathway of the yeast cell. It decreases the concentration of acetyl-CoA significantly. On the other hand, if we set the glucose inhibition as going to zero and the inhibition by acetaldehyde as high enough then we get the simulation results as shown in Figure 3. We can observe that these conditions influence not only the concentration of pyruvate but also the concentration of the other metabolites including the concentration of ethanol. The maximum concentration of ethanol decreases as the acetaldehyde inhibition increases. However, when the acetaldehyde inhibition is high enough, it did not influence the concentration of acetyl-CoA as the main substrate for cell respiration (see Figure 3 (C)). It means that the feedback inhibitions by glucose significantly affect the respiratory pathway of the yeast cell, whereas acetaldehyde inhibition significantly affects the fermentation pathway. This is an interesting result since the two pathways, respiration and fermentation pathways, are important pathways for yeast cells that must be maintained.

Sensitivity analysis using metabolic control analysis
In this section we analyse the sensitivity of reactions involved in the fermentation pathway. Metabolic control analysis is a quantitative tool that quantifies the effect of perturbation in a metabolic network. It also quantifies the importance of each enzyme in controlling the flux of ethanol. In the metabolic control analysis, the sensitivity of reactions is measured via the following coefficients, namely elasticity coefficient and control coefficient for flux and metabolite. The elasticity coefficient quantifies the effect of small perturbation of a reaction parameter on the local reaction, whereas the control coefficient quantifies global effects of the perturbation of enzyme activity to the flux or metabolite concentration. To derive the sensitivity analysis of our system, now consider the nonlinear kinetic model in Equation (7). In a simple form, the kinetic model (7) can be rewritten in the following form: where r ¼ r 1 ; . . . ; r 10 ð Þ ; S ¼ ðS 1 ; . . . ; S 6 Þ; p is a vector of kinetic parameter of model (7), and N is a 6 Â 10-dimensional matrix of stoichiometric coefficients (the stoichiometric matrix) of system (7), Table 1. The kinetic parameters of system (7) taken from [7].
The row echelon form of the augmented matrix NjI ð Þ is given as the following form: (25) (26) r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 9 r 10 _ It can be identified that there are no zero rows in its row echelon form meaning that no reduction is necessary for metabolites of system (24) since there is no dependence of the differential equations. Now let us investigate the flux of system (24). The flux of system (24) is defined as the steady-state reaction rate that can be reached when Nr S; p ð Þ ¼ 0: Therefore, finding the flux of the system (24) equals with finding the kernel of the Stoichiometric matrix N. Let J i denote the flux of reaction r i : Since the rank of the matrix N is 6, there exist four independent fluxes that can be chosen freely. Choosing J 4 ; . . . ; J 7 as the independent  fluxes (the r 4 , … ,r 7 -columns without pivots in the echelon form of the stoichiometric matrix N), we then have with the flux vector J ¼ ðJ 1 ; . . . ; J 10 Þ; the independent flux vector J id ¼ ðJ 4 ; . . . ; J 7 Þ and the kernel of matrix N, highest negative control is given by reaction r 5 catalysed by acetaldehyde dehydrogenase. These results are in agreement with our previous analysis results [9]. It means that, if we increase the supply of glucose, then it will affect the increasing of ethanol production. However, the increasing of glucose supply should be under control since high glucose concentration significantly affects the respiratory pathway of the yeast cell. Otherwise, if we decrease reaction r 5 catalysed by acetaldehyde dehydrogenase then it also will increase the production of ethanol. Some simulations are presented to confirm this metabolic control analysis results. For instance, in Figure 5, we present the comparison of several regulations according to the sensitivity results. We can observe that when the glucose supply is increased about 100% (reaction with positive effect), it also increases the maximum concentration of ethanol for its steady state concentration. However, the increasing of ethanol concentration is not really significant. In comparison with the increasing of the maximum reaction rate of r 4 by about 50% (reaction with positive effect), this regulation increases the concentration of ethanol higher than the previous regulation. Moreover, when the regulation is focused on the reaction r 5 by decreasing its maximal reaction rate  (reaction with negative effect), then this regulation significantly affects the increasing of ethanol concentration compared with the previous two regulations. The maximal concentration of ethanol increases dramatically, i.e. almost more than 100%.
If we combine several regulation rules to increase the production of ethanol, we get the simulation results depicted in Figure 6. When the maximal reaction rates of r 4 and r 5 are simultaneously regulated we find that the maximal concentration of ethanol increases twice from the previous regulation (regulating only r 4 ). When the maximal reaction rates of r 4 , r 5 and glucose supply are simultaneously regulated we find that the maximum concentration of ethanol increases slightly compared with regulation of r 4 and r 5 . Since high glucose concentration will affect the inhibition at r 3 and r 6 then the best regulation choice is by regulating r 4 and r 5 simultaneously.

Conclusions
The kinetic model of the main metabolism of fermentation system was modelled by taking into consideration feedback inhibition effects of glucose and acetaldehyde to some reactions. The inhibition mechanisms occurred competitively and noncompetitively. The model produced at least one stable steady state and at most two stable steady states, if it fulfilled certain conditions. Numerical observations showed that high inhibition by glucose influenced the concentration of acetyl-CoA that will lead to the respiration pathway of the yeast cell, whereas high inhibition of acetaldehyde affected the concentration of all metabolites including the maximal concentration of ethanol. By applying metabolic control analysis to the metabolic system we found that there is one regulation rule that can be applied to increase the production of ethanol. Based on our numerical observations we found that the best regulation rule was by regulating the two key reactions that were catalysed by alcohol dehydrogenase and acetaldehyde dehydrogenase. By regulating the maximal activity of the two enzymes simultaneously, the maximal concentration of ethanol can be increased significantly.