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.

Keywords:

• Yeast cell
• fermentation system
• feedback inhibition
• kinetic modelling
• stability theory
• metabolic control analysis
• metabolite dynamics

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–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–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.

Methods

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₁) is converted into pyruvate (S₂) 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(7a) $$\dot{S_1}\left(t\right)=\mathrm{ }\alpha G-\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1},$$(7a) ,9](9) $$S_1^{\mathrm{*}}=\frac{\alpha G{K}_1}{V_1-\alpha G\left(K_a{S}_4^{\mathrm{*}}+1\right)}.$$(9) . Since there exists uncompetitive internal inhibition by acetaldehyde (S₄) to the glycolytic pathway [20], the Michaelis-Menten kinetic equation becomes (1) $$r_1=V_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1},$$(1) where $K_{a }$is an inhibition constant of $S_4$with respect to $r_1,$$\mathrm{ }{V}_1$ is the maximum reaction of glucose conversion to become pyruvate in the glycolysis pathway and K₁ 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 (2) $$r_2=V_2\frac{S_2}{S_2+K_2},$$(2) 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: (3) $$r_3=V_3\frac{S_2}{S_2+K_3}·\frac{1}{K_b{S}_1+1},$$(3) 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₅) and ethanol (S₆). 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: (4) $$r_4=\mathrm{ }\frac{{V_{4f}{K}_{4b}S}_4-{V_{4b}{K}_{4f}S}_6}{{K_{4f}K}_{4b}+{K_{4b}S}_4+{K_{4f}S}_6},$$(4) (5) $${\mathrm{r}}_5={\mathrm{V}}_5\frac{{\mathrm{S}}_4}{{\mathrm{S}}_4+{\mathrm{K}}_5}.$$(5)

Figure 1. Central metabolism of Saccharomyces cerevisiae. Reaction r₁ is a lumped process of glucose conversion to become pyruvate and reaction r_i, i = 2,…,6 is a chemical conversion catalysed by: (2) pyruvate decarboxylase, (3) pyruvate dehydrogenase complex, (4) alcohol dehydrogenase, (5) acetaldehyde dehydrogenase and (6) acetyl-CoA synthetase.

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: (6) $$r_6=V_6\frac{S_5}{S_5+K_6}·\frac{1}{K_c{S}_1+1},$$(6) 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₆ 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: (7a) $$\dot{S_1}\left(t\right)=\mathrm{ }\alpha G-\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1},$$(7a) (7b) $$\dot{S_2}\left(t\right)=\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1}-\left(V_2\frac{S_2}{S_2+K_2}+V_3\frac{S_2}{S_2+K_3}·\frac{1}{K_{\mathrm{b}}{S}_1+1}\right),$$(7b) (7c) $$\dot{S_3}\left(t\right)=\mathrm{ }{V}_3\frac{S_2}{S_2+K_3}·\frac{1}{K_{\mathrm{b}}{S}_1+1}+V_6\frac{S_5}{S_5+K_6}·\frac{1}{K_{\mathrm{c}}{S}_1+1}-\mathrm{ }\alpha S_3,$$(7c) (7d) $$\dot{S_4}\left(t\right)=\mathrm{ }{V}_2\frac{S_2}{S_2+K_2}-\mathrm{ }\left(\frac{{V_{4f}{K}_{4b}S}_4-{V_{4b}{K}_{4f}S}_6}{{K_{4f}K}_{4b}+{K_{4b}S}_4+{K_{4f}S}_6}+\mathrm{ }{V}_5\frac{S_4}{S_4+K_5}\right),$$(7d) (7e) $$\dot{S_5}\left(t\right)=\mathrm{ }{V}_5\frac{S_4}{S_4+K_5}-V_6\frac{S_5}{S_5+K_6}·\frac{1}{K_{\mathrm{c}}{S}_1+1}-\mathrm{ }\alpha S_5,$$(7e) (7f) $$\dot{S_6}\left(t\right)=\mathrm{ }\frac{{V_{4f}{K}_{4b}S}_4-{V_{4b}{K}_{4f}S}_6}{{K_{4f}K}_{4b}+{K_{4b}S}_4+{K_{4f}S}_6}-\mathrm{ }\alpha S_6,$$(7f) with initial conditions (7g) $$S_1\left(0\right)=S_a{,\mathrm{ }S}_2\left(0\right)=S_b,\mathrm{ }{S}_3\left(0\right)=S_c,S_4\left(0\right)=S_d,S_5\left(0\right)=S_e,S_6\left(0\right)=S_f,$$(7g) where $S_a,$$\dots ,$ $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 ${\mathbb{R}}_{+}^6=\left\{\left(S_1,\dots ,S_6\right)\in {\mathbb{R}}^6|S_i\ge 0,\mathrm{ }i=1,\dots ,6\right\}.$

At the equilibrium condition Equation (7(7a) $$\dot{S_1}\left(t\right)=\mathrm{ }\alpha G-\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1},$$(7a) a) becomes (8) $$\alpha G-\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1}=0$$(8)

such that we get the equilibrium solution for S₁ depending on the equilibrium solution of S₄, (9) $$S_1^{\mathrm{*}}=\frac{\alpha G{K}_1}{V_1-\alpha G\left(K_a{S}_4^{\mathrm{*}}+1\right)}.$$(9)

The equilibrium solution of S₁ will be positive if it fulfils the conditions (10) $$V_1>\alpha \mathrm{G\hspace{0.5em}and\hspace{0.5em}}0<S_4^{\mathrm{*}}<\frac{V_1\mathrm{ }-\alpha G}{\alpha G{K}_a}.$$(10)

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₄ should be maintained in the interval $\left(0,\frac{V_1-\alpha G}{\alpha G{K}_a}\right).$ From Equation (7(7a) $$\dot{S_1}\left(t\right)=\mathrm{ }\alpha G-\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1},$$(7a) b), we get (11) $$V_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1}-\left(V_2\frac{S_2}{S_2+K_2}+V_3\frac{S_2}{S_2+K_3}·\frac{1}{K_b{S}_1+1}\right)=0,$$(11) which can be rewritten in quadratic polynomial, (12) $$f\left(S_2\right)=\mathrm{ }{a}_1{S_2}^2+a_2{S}_2+\mathrm{ }{a}_3=0,$$(12) with $$a_1=\mathrm{ }{V}_1{S}_1{(K}_b{S}_1+1)-\left[S_1\left(K_a{S}_4+1\right)+K_1\right]\left[{V_2(K}_b{S}_1+1\right)+V_3],$$ $$a_2=\mathrm{ }{V_1{S}_1{(K}_b{S}_1+1\left)\right(K}_2+K_3)-\left[S_1\left(K_a{S}_4+1\right)+K_1\right]{[V}_2{K}_3{(K}_b{S}_1+1)+V_3{K}_2],$$ $$a_3={V_1{S}_1{(K}_b{S}_1+1)K}_2{K}_3.$$

The roots of polynomial (12) are equilibrium solutions of S₂ which depend on the equilibrium solutions of S₁ and S₄. If the equilibrium solution of S₁ fulfils condition (10), then we have a₃>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₂. 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₃: (13) $$S_3^{\mathrm{*}}\mathrm{ }=\frac{V_3{S}_2\left(S_5+K_6\right)\left(K_{\mathrm{c}}{S}_1+1\right)\mathrm{ }+V_6{S}_5\left(S_2+K_3\right)\left(K_{\mathrm{b}}{S}_1+1\right)}{\alpha \left(S_2+K_3\right)\left(K_{\mathrm{b}}{S}_1+1\right)\left(S_5+K_6\right)\left(K_{\mathrm{c}}{S}_1+1\right)}$$(13)

that will be positive if the equilibrium solutions of S₁, S₂ and S₅ are positive. From (7d), we have (14) $$V_2\frac{S_2}{S_2+K_2}-\mathrm{ }\left(\frac{{V_{4f}{K}_{4b}S}_4-{V_{4b}{K}_{4f}S}_6}{{K_{4f}K}_{4b}+{K_{4b}S}_4+{K_{4f}S}_6}+\mathrm{ }{V}_5\frac{S_4}{S_4+K_5}\right)=0,$$(14)

that can be rewritten as a quadratic polynomial as follows: (15) $$f\left(S_4\right)=\mathrm{ }{b}_1{S_4}^2+b_2{S}_4+\mathrm{ }{b}_3=0,$$(15) with $$b_1=V_2{S}_2{K}_{4b}\mathrm{ }–\mathrm{ }\left[{\left(S_2+K_2\right)V}_{4f}{K}_{4b}+\mathrm{ }{\left(S_2+K_2\right)K}_{4b}{V}_5\right],$$ $$b_2=\mathrm{ }{V}_2{S}_2\left({K_{4f}K}_{4b}+{K_{4f}S}_6\right)+\mathrm{ }{V_2{S}_{2}K}_{4b}{K}_5+{{\left(S_2+K_2\right)V}_{4b}{K}_{4f}S}_6-\mathrm{ }\left[{\left(S_2+K_2\right)V}_{4f}{K}_{4b}{K}_5+\left(S_2+K_2\right){K_{4f}S}_6{V}_5+{\left(S_2+K_2\right)V}_{4f}{K}_{4b}{K}_5\right],$$ $$b_3=\mathrm{ }{V}_2{S}_2{K}_5\left({K_{4f}K}_{4b}+{K_{4f}S}_6\right)+{{\left(S_2+K_2\right)V}_{4b}{K}_{4f}S}_6{K}_5.$$

If the equilibrium solutions of S₂ and S₆ are positive, then we get b₃>0. Using the similar analysis with S₂, we find that the polynomial (15) will produce at most two positive roots for S₄. This result also leads to the possibility of multiple steady states in the modelled fermentation system. Furthermore, from (7e) we have (16) $$V_5\frac{S_4}{S_4+K_5}-V_6\frac{S_5}{S_5+K_6}·\frac{1}{K_{\mathrm{c}}{S}_1+1}-\mathrm{ }\alpha S_5=0,$$(16)

that also can be rewritten in a quadratic polynomial form, i.e. (17) $$f\left(S_5\right)=\mathrm{ }{c}_1{S_5}^2+c_2{S}_5+\mathrm{ }{c}_3=0,$$(17) with $$c_1=\alpha \left(S_4+K_5\right)\left(K_c{S}_1+1\right),$$ $$c_2= \left[\alpha K_6\left(S_4+K_5\right)\left(K_c{S}_1+1\right)\right]-V_5{S}_4\left(K_c{S}_1+1\right),$$ $$c_3= -V_5{S}_4{K}_6\left(K_c{S}_1+1\right).$$

Since c₁>0 and c₃<0 when S₁ and S₄ are positive, then polynomial (17) will produce at most one positive equilibrium solution for S₅.

Furthermore, from the last Equation (7(7a) $$\dot{S_1}\left(t\right)=\mathrm{ }\alpha G-\mathrm{ }{V}_1\frac{S_1}{S_1\left(K_a{S}_4+1\right)+K_1},$$(7a) f), we have (18) $$\frac{{V_{4f}{K}_{4b}S}_4-{V_{4b}{K}_{4f}S}_6}{{K_{4f}K}_{4b}+{K_{4b}S}_4+{K_{4f}S}_6}-\mathrm{ }\alpha S_6=0.$$(18)

In quadratic polynomial form, Equation (18)(18) $$\frac{{V_{4f}{K}_{4b}S}_4-{V_{4b}{K}_{4f}S}_6}{{K_{4f}K}_{4b}+{K_{4b}S}_4+{K_{4f}S}_6}-\mathrm{ }\alpha S_6=0.$$(18) can be rewritten as follows: (19) $$f\left(S_6\right)=\mathrm{ }{d}_1{S_6}^2+d_2{S}_6+\mathrm{ }{d}_3=0,$$(19) with $${d_1=K}_{4f}\alpha ,$$ $$d_2={K_{4f}K}_{4b}\alpha +V_{4b}{K}_{4f}+{K_{4b}S}_4\alpha ,$$ $$d_3=-{V_{4f}{K}_{4b}S}_4.$$

Since d₁>0 and d₃<0 when S₄>0, we get that S₄ has only one positive equilibrium solution.

Summarising our analysis results, we found that the fermentation system (7) has an equilibrium solution $E=\left({S_1}^{\mathrm{*}},\mathrm{ }{S_{2i}}^{\mathrm{*}},\mathrm{ }{S_3}^{\mathrm{*}},{S_{4i}}^{\mathrm{*}},{S_5}^{\mathrm{*}},{S_6}^{\mathrm{*}}\right)$ where ${S_{2i}}^{\mathrm{*}}$and ${S_{4i}}^{\mathrm{*}},$ respectively, stand for the i-th positive root of S₂ and S₄. 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: (20) $$J=\left(\begin{array}{cc}\begin{array}{ccc}{-\Delta }_1{K}_1& 0& 0\\ {\Delta }_{10}& -{\Delta }_9& 0\\ -{\Delta }_8& {\Delta }_2{\Delta }_7{K}_3& -\alpha \end{array}& \begin{array}{ccc}{\Delta }_1{{\left({S_1}^{\mathrm{*}}\right)}^{2}K}_a& 0& 0\\ -{\Delta }_1{\left({S_1}^{\mathrm{*}}\right)}^2{K}_a& 0& 0\\ 0& \frac{{\Delta }_3{K}_6}{{S_5}^{\mathrm{*}}+K_6}& 0\end{array}\\ \begin{array}{ccc}0& {{\Delta }_6}^2{V}_2{K}_2& 0\\ {\Delta }_{12}& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}-{\Delta }_{11}& \mathrm{ }0& {\Delta }_5\\ \frac{V_5{K}_5}{{\left({S_4}^{\mathrm{*}}+K_5\right)}^2}& -(\frac{{\Delta }_3{K}_6}{{S_5}^{\mathrm{*}}+K_6}+\alpha )& 0\\ {\Delta }_4& 0& -({\Delta }_5+\alpha )\end{array}\end{array}\right),$$(20) where $${\Delta }_1=\frac{V_1}{{\left[{S_1}^{\mathrm{*}}\left(K_a{S_4}^{\mathrm{*}}+1\right)+K_1\right]}^2},$$ $${\Delta }_2=\frac{V_3}{\left({S_2}^{\mathrm{*}}+K_3\right)\left(K_b{S_1}^{\mathrm{*}}+1\right)},$$ $${\Delta }_3=\frac{V_6}{\left({S_5}^{\mathrm{*}}+K_6\right)\left(K_c{S_1}^{\mathrm{*}}+1\right)},$$ $${\Delta }_4=\frac{K_{4f}{K}_{4b}\left(V_{4f}{K}_{4b}+V_{4f}{S_6}^{\mathrm{*}}+V_{4b}{S_6}^{\mathrm{*}}\right)}{{\left(K_{4f}{K}_{4b}+K_{4b}{S_4}^{\mathrm{*}}+K_{4f}{S_6}^{\mathrm{*}}\right)}^2},$$ $${\Delta }_5=\frac{K_{4f}{K}_{4b}\left(V_{4b}{K}_{4f}+V_{4b}{S_4}^{\mathrm{*}}+V_{4f}{S_4}^{\mathrm{*}}\right)}{{\left(K_{4f}{K}_{4b}+K_{4b}{S_4}^{\mathrm{*}}+K_{4f}{S_6}^{\mathrm{*}}\right)}^2},$$ $${\Delta }_6=\frac{1}{{S_2}^{\mathrm{*}}+K_2},$$ $${\Delta }_7=\frac{1}{{S_2}^{\mathrm{*}}+K_3},$$ $${\Delta }_8=\frac{K_b{S_2}^{\mathrm{*}}{\Delta }_2}{K_b{S_1}^{\mathrm{*}}+1}+\frac{{S_5}^{\mathrm{*}}{K}_c{\Delta }_3}{1+K_c{S_1}^{\mathrm{*}}},$$ $${{\Delta }_9=V}_2{{\Delta }_6}^2{K}_2+{\Delta }_7{\Delta }_2{K}_3,$$ $${{\Delta }_{10}=\Delta }_1{K}_1+{\Delta }_2{\left({S_2}^{\mathrm{*}}\right)}^2{K}_b,$$ $${{\Delta }_{11}=A}_4+\frac{V_5{K}_5}{{\left({S_4}^{\mathrm{*}}+K_5\right)}^2},$$ $${\Delta }_{12}=\frac{{\Delta }_3{S_5}^{\mathrm{*}}{K}_c}{1+K_c{S_1}^{\mathrm{*}}}.$$

We can observe that two eigenvalues of matrix J are negative, i.e. (21) $${\lambda }_1=-\left(\frac{{\Delta }_3{K}_6}{S_5+K_6}+\alpha \right){\text{\hspace{0.5em}and\hspace{0.5em}}\lambda }_2=-\alpha .$$(21)

The rest of eigenvalues of matrix J is the solutions of the characteristic equation, (22) $${\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0,$$(22) where coefficients $a_1,\dots , a_4$ depend on ${\mathrm{\Delta }}_1,\dots ,{\mathrm{\Delta }}_{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

1. $a_1>0,$

2. $\left|\begin{array}{cc}{a}_1& 1\\ 0& a_2\end{array}\right|>0$ or $a_1{a}_2>0,$

3. $\left|\begin{array}{ccc}{a}_1& 1& 0\\ a_3& a_2& a_1\\ 0& 0& a_3\end{array}\right|>0$ or $a_1{a}_2{a}_3-{a_3}^2>0,$

4. $\left|\begin{array}{cc}\begin{array}{cc}{a}_1& 1\\ a_3& a_2\end{array}& \begin{array}{cc}0& 0\\ a_1& 1\end{array}\\ \begin{array}{cc}0& a_4\\ 0& 0\end{array}& \begin{array}{cc}{a}_3& a_2\\ a_4& a_3\end{array}\end{array}\right|>0$ or ${a_3}^3+a_1{a}_2{a_3}^2-{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.

Table 1. The kinetic parameters of system (7) taken from [7].

Par.	Val.	Unit	Par.	Val.	Unit	Par.	Val.	Unit
$v_1$	1, 43	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{1\mathrm{ }}$	0, 94	$\mathrm{g}\mathrm{ }{\mathrm{L}}^{-1}$	$\alpha$	0,01	${\mathrm{h}}^{-1}$
$v_2$	5, 81	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{2\mathrm{ }}$	5 ×${10}^{-7}$	$\mathrm{g}{\text{ L}}^{-1}$	G	15	$\mathrm{g}\mathrm{ }{\mathrm{L}}^{-1}$
$v_3$	0, 501	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{3\mathrm{ }}$	2 ×${10}^{-5}$	$\mathrm{g}{\text{ L}}^{-1}$	$k_{1i}$	14, 2	$\mathrm{g}{\text{ L}}^{-1}$
$v_{4f}$	2, 82	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{4f}$	0,034	$\mathrm{g}{\text{ L}}^{-1}$	$k_{3i}$	0, 101	$\mathrm{g}{\text{ L}}^{-1}$
$v_{4b}$	0,0125	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{4b}$	0,057	$\mathrm{g}{\text{ L}}^{-1}$	$k_{6i}$	440	$\mathrm{g}{\text{ L}}^{-1}$
$v_5$	4,8	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{5\mathrm{ }}$	2,64 ×${10}^{-4}$	$\mathrm{g}{\text{ L}}^{-1}$
$v_6$	0, 0104	$\mathrm{g}{\mathrm{g}}^{-1}{\mathrm{h}}^{-1}$	$k_{6\mathrm{ }}$	0, 0102	$\mathrm{g}{\text{ L}}^{-1}$

Results and discussion

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) $$E=(S_1^{\mathrm{*}},S_2^{\mathrm{*}},S_3^{\mathrm{*}},S_4^{\mathrm{*}},S_5^{\mathrm{*}},S_6^{\mathrm{*}})\text{=}(0.1102;\text{1.3}{\times 10}^{-8};\; 0.0439;\; 8.4\times {10}^{-6};\; 14.9432;\; 0.0031)$$(23) with eigenvalues ${\mathrm{\lambda }}_1=-0.01,$ ${\mathrm{\lambda }}_2=-0.01,$ ${\mathrm{\lambda }}_3=-1.22, {\lambda }_4=-1.103{\times 10}^7,{\lambda }_5=-17142.3, {\lambda }_6=-0.22$ implying the local stability of E. By choosing initial conditions S₁(0)=20 g L⁻¹ and zero for the others and final time observation 100 hours, we get numerical solutions of system (7) depicted in Figure 2.

Figure 2. Comparison of the transient behaviour of metabolites in the central metabolism of yeast cells with and without inhibition effects.

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.

Figure 3. Transient behaviour of metabolites when the glucose inhibition is over and the inhibition by acetaldehyde is very high.

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: (24) $$\frac{d\mathit{S}}{\mathit{dt}}=N\mathit{r}\left(\mathit{S},\mathit{p}\right),$$(24) where $\mathit{r}=\left(r_1, \dots , r_{10}\right),\mathit{S}=(S_1,\dots ,S_6),\mathit{p}$ is a vector of kinetic parameter of model (7), and N is a $6\times 10$-dimensional matrix of stoichiometric coefficients (the stoichiometric matrix) of system (7), (25) $$N=\left(\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}-1\\ 1\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ -1\\ 0\end{array}\\ \begin{array}{c}1\\ 0\\ 0\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ -1\\ 1\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}-1\\ 0\\ 1\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}-1\\ 1\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\\ 1\end{array}\\ \begin{array}{c}0\\ -1\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ 0\\ -1\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ -1\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ -1\end{array}\end{array}\end{array}\end{array}\end{array}\right).$$(25)

Reactions $r_1, \dots , r_6$ are given in Equations (1)–(6), whereas reactions $r_7, \dots , r_{10},$respectively, represent the feeding rate of glucose and the outflows of the metabolite products, i.e. $r_7=\alpha \mathrm{G},$ $r_8=\alpha S_3,r_9=\alpha S_5,r_{10}=\alpha S_6.$ If we extend the stoichiometric matrix N using an $6\times 6$-dimensional identity matrix then we get the augmented matrix$\left(N\right|I)$where each column of the identity matrix I represents a time derivative of system (24), i.e. (26) $$\left(N|I\right)=\left(\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}-1\\ 1\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ -1\\ 0\end{array}\\ \begin{array}{c}1\\ 0\\ 0\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ -1\\ 1\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}-1\\ 0\\ 1\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}-1\\ 1\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\\ 1\end{array}\\ \begin{array}{c}0\\ -1\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ 0\\ -1\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ -1\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ -1\end{array}\end{array}\end{array}\end{array}\end{array}|\begin{array}{cc}\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}& \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\end{array}\right).$$(26)

The row echelon form of the augmented matrix $\left(N|I\right)$ is given as the following form:

Table

	$r_1$	$r_2$	$r_3$	$r_4$	$r_5$	$r_6$	$r_7$	$r_8$	$r_9$	$r_{10}$	$\dot{S_1}$	$\dot{S_2}$	$\dot{S_3}$	$\dot{S_4}$	$\dot{S_5}$	$\dot{S_6}$
$S_1$	–1	0	0	0	0	1	0	0	0	0	1	0	0	0	0	0
$S_2$	0	0	–1	–1	–1	0	1	0	0	0	1	1	0	1	0	0
$S_3$	0	0	0	−1	–1	1	1	–1	0	0	1	1	1	1	0	0
$S_4$	0	1	0	1	–1	0	0	0	0	0	0	0	0	1	0	0
$S_5$	0	0	0	0	1	–1	0	0	–1	0	0	0	0	0	1	0
$S_6$	0	0	0	1	0	0	0	0	0	–1	0	0	0	0	0	1

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 $N\mathit{r}\left(\mathit{S},\mathit{p}\right)=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,\dots , J_7$ as the independent fluxes (the r₄,…,r₇-columns without pivots in the echelon form of the stoichiometric matrix N), we then have (27) $$\mathit{J}=\mathit{K}{\mathit{J}}_{\mathit{id}},$$(27) with the flux vector $\mathit{J}=(J_1,\dots , J_{10}),$ the independent flux vector ${\mathit{J}}_{\mathit{id}}=(J_4,\dots , J_7)$ and the kernel of matrix N, (28) $$\mathit{K}=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 1\\ \begin{array}{c}-1\\ 1\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ -1\\ \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ 1\\ \begin{array}{c}-1\\ 0\end{array}\end{array}\\ \begin{array}{c}1\\ 0\\ \begin{array}{c}0\\ -1\\ \begin{array}{c}1\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ 1\\ \begin{array}{c}0\\ 1\\ \begin{array}{c}-1\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ \begin{array}{c}1\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}1\\ 1\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right].$$(28)

The matrix K shows the dependence of the steady state fluxes on the independent fluxes. In a normal form we have (29) $$\mathit{Ҡ}= {D_{\mathit{J}}}^{-1}\mathit{K}{D}_{{\mathit{J}}_{\mathit{id}}}=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ \frac{J_4}{J_2}\\ \begin{array}{c}-\frac{J_4}{J_3}\\ \frac{J_4}{J_4}\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ -\frac{J_4}{J_8}\\ \begin{array}{c}0\\ -\frac{J_4}{J_{10}}\end{array}\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ \frac{J_5}{J_2}\\ \begin{array}{c}-\frac{J_5}{J_3}\\ 0\end{array}\end{array}\\ \begin{array}{c}\frac{J_5}{J_5}\\ 0\\ \begin{array}{c}0\\ -\frac{J_5}{J_8}\\ \begin{array}{c}\frac{J_5}{J_9}\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ \frac{J_6}{J_6}\\ \begin{array}{c}0\\ \frac{J_6}{J_8}\\ \begin{array}{c}-\frac{J_6}{J_9}\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\frac{J_7}{J_1}\\ \frac{-2J_7}{J_2}\\ \begin{array}{c}\frac{J_7}{J_3}\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}\frac{J_7}{J_7}\\ \frac{J_7}{J_8}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right],$$(29) where $D_{\mathit{J}}$ and $D_{{\mathit{J}}_{\mathit{id}}}$ stand for the diagonal matrix of flux vector J and independent flux ${\mathit{J}}_{\mathit{id}},$ respectively.

In metabolic control theory, there exist two theorems that can be applied to find the matrix control for flux and metabolite of the system (for detail about the two theorems see [23]). The normalised form of both theorems is written as (30) $$\left(\begin{array}{c}{C}_J\\ C_S\end{array}\right)\left(\begin{array}{cc}\mathrm{Ҡ}& {-M}_{\epsilon }\end{array}\right)=\left(\begin{array}{cc}\mathrm{Ҡ}& 0\\ 0& I\end{array}\right),$$(30) where $C_J$ and $C_S,$ respectively, represent the normalised matrix control for flux and metabolite, $\mathrm{Ҡ}$ is the normalised form of the kernel of N given in Equation (29)(29) $$\mathit{Ҡ}= {D_{\mathit{J}}}^{-1}\mathit{K}{D}_{{\mathit{J}}_{\mathit{id}}}=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ \frac{J_4}{J_2}\\ \begin{array}{c}-\frac{J_4}{J_3}\\ \frac{J_4}{J_4}\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ -\frac{J_4}{J_8}\\ \begin{array}{c}0\\ -\frac{J_4}{J_{10}}\end{array}\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ \frac{J_5}{J_2}\\ \begin{array}{c}-\frac{J_5}{J_3}\\ 0\end{array}\end{array}\\ \begin{array}{c}\frac{J_5}{J_5}\\ 0\\ \begin{array}{c}0\\ -\frac{J_5}{J_8}\\ \begin{array}{c}\frac{J_5}{J_9}\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ \frac{J_6}{J_6}\\ \begin{array}{c}0\\ \frac{J_6}{J_8}\\ \begin{array}{c}-\frac{J_6}{J_9}\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\frac{J_7}{J_1}\\ \frac{-2J_7}{J_2}\\ \begin{array}{c}\frac{J_7}{J_3}\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}\frac{J_7}{J_7}\\ \frac{J_7}{J_8}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right],$$(29) , I is an identity matrix, and $M_{\epsilon }$ is a matrix whose entries are elasticity coefficients, ${\epsilon }_s^r=\frac{s}{r}\frac{\partial r}{\partial s}.$ By differentiating reaction rate $r_i, i=1,.,10$ with respect to the metabolite variable $S_i, i=1,\dots ,6,$ we get the elasticity matrix as follows: (31) $$M_{\epsilon }=\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}\frac{K_1}{m_1}\\ 0\\ \begin{array}{c}\frac{-S_1{K}_{3i}}{\left(S_1{K}_{3i}+1\right)}\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ \frac{-S_1{K}_{6i}}{\left(S_1{K}_{6i}+1\right)}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ \frac{K_2}{S_2+K_2}\\ \begin{array}{c}\frac{K_3}{S_2+K_3}\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\\ \begin{array}{c}0\\ 1\\ \begin{array}{c}0\\ 1\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}\frac{-S_1{K}_{1i}}{m_1}\\ 0\\ \begin{array}{c}0\\ \frac{m_2}{m_3}\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\\ \begin{array}{c}\frac{K_5}{S_5+K_5}\\ \frac{K_6}{S_5+K_6}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}1\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ \frac{m_4}{m_3}\end{array}\end{array}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right].$$(31)

Therefore, by substituting matrix $Ҡ$ in (29) and matrix $M_{\epsilon }$ in (31) into Equation (30)(30) $$\left(\begin{array}{c}{C}_J\\ C_S\end{array}\right)\left(\begin{array}{cc}\mathrm{Ҡ}& {-M}_{\epsilon }\end{array}\right)=\left(\begin{array}{cc}\mathrm{Ҡ}& 0\\ 0& I\end{array}\right),$$(30) , we can derive the matrix control for flux, $C_J,$ and the matrix control for metabolites, $C_S.$

Using the kinetic parameters in Table 1, the stable positive steady state E, and Equations (29)–(31), we get the values of the sensitivity of ethanol flux ($J_{10})$ in term of flux control coefficient ($C_{r_k}^{J_{10}}$) which are depicted in Figure 4. Here, we only present the sensitivity of ethanol flux ($J_{10}$) with respect to the reaction rate changes, $r_k, k=1,\dots ,10,$ because our analysis is focused on the effects of the existence of feedback inhibition in some reactions to the ethanol production.

Figure 4. Distribution of flux control coefficients for quantifying the sensitivity of ethanol flux with respect to the reaction rate changes.

From the metabolic control analysis results, we can observe that each reaction has a different effect to the ethanol flux. Some reactions have almost zero effect to the ethanol flux including reactions that experience feedback inhibitions, whereas other reactions have positive effect and negative effect. The highest positive control is given by the glucose supply with control value $C_{r_7}^{J_{10}}=1,\mathrm{ }0296,$ whereas the highest negative control is given by reaction r₅ 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₅ 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₄ 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₅ 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%.

Figure 5. Comparison of ethanol concentration for different regulation rules based on the metabolic control analysis results.

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₄ and r₅ are simultaneously regulated we find that the maximal concentration of ethanol increases twice from the previous regulation (regulating only r₄). When the maximal reaction rates of r₄, r₅ and glucose supply are simultaneously regulated we find that the maximum concentration of ethanol increases slightly compared with regulation of r₄ and r₅. Since high glucose concentration will affect the inhibition at r₃ and r₆ then the best regulation choice is by regulating r₄ and r₅ simultaneously.

Figure 6. Comparison of ethanol concentration for different combination of regulation rules based on the metabolic control analysis results.

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.

Additional information

Funding

The authors would like to thank DIKTI via Institute for Research and Community Services Hasanuddin University that supported part of this research through Penelitian Dasar Unggulan Perguruan Tinggi (PDUPT) under research grant No. 2610/UN4.21/LK.23/2017 and 1740/UN4.21/PL.01.10/2019.