Abstract: In recent years constraint based methods such as Flux Balance Analysis (FBA) has widely applied for computation of flux distributions in the metabolic networks. In this study, the effects of changes in the intercellular concentration on the Gibbs free energy of the system and subsequently on the model’s fluxes have studied. This method which makes a correlation between the flux directions and metabolite concentrations, has applied to large scale metabolic network, Escherichia coli iAF1260. The biomass and succinate fluxes have selected as objective functions and the multi objective genetic algorithm has used to optimization of the E. coli iAF1260 network. The obtained results revealed that a living system such as E. coli and its mutants are not stable from thermodynamic point of view before reaching to the cellular threshold concentrations. Also the behavior of the mutants of microorganism involving their return to the wild type phenotypes could be justified.
INTRODUCTION
Reconstructed metabolic networks in living organisms are going to be a powerful tool in the prediction of the phenotypes and analyzing their sensitivity to the imposed environmental changes and internal stimulus. Genome-scale metabolic networks which have reconstructed for some microorganisms such as Escherichia coli and Saccharomyces cerevisiae and their high abilities for predicting several phenotypes (in agreement with experimental data) have previously demonstrated (Lee et al., 2002; Salgado et al., 2004; Shen-Orr et al., 2002; Guelzim et al., 2002). For example, Genome-scale metabolic networks in Escherichia coli and saccharomyces cerevisiae have predicted 86 and 88% of their genes deletion phenotypes, respectively (Natalie et al., 2004). Flux Balance Analysis (FBA) is one of the most important tools in order to study of the metabolic networks. It should be noted that the FBA is a fair steady state approximation for estimating the behavior of the cells when there is not enough kinetic information about the reactions of the metabolic networks (Varma et al., 1993; Varma and Palsson, 1993; Edwards et al., 2001; Dien and Lidstrom, 2002; Kayser et al., 2005). Genome-scale metabolic networks can be reconstructed through the definition of the stoichiometry of matrix from mass balance equations and then imposing biological and systemic constraints into the network, determining the objective functions, energy requirements and other cell demands (Covert et al., 2004; Price et al., 2004). In this regard, the constraint-based analysis is a conventional method to analyse the metabolic networks and provide a framework to compute the cellular functions (Palsson, 2000). In this method that based on the convex analysis by sequential imposing the biological and systemic constraints in the network, the allowable solution space of phenotypes is gradually shrinks. This modeling method considers cellular, biochemical and systemic functions limits in the allowable solution space.
The important point in the of use of the constraint base analysis is that the results of such analysis let us know what the networks could do rather than determining the exact network phenotypes with respect to the specific external environment. Imposing additional constraints (physiologic, metabolic and thermodynamic constraints) in FBA restrict the number of possible solutions. These solutions are more compatible with real accessible solutions in metabolic networks.
Also, the Energy Balance Analysis (EBA) which incorporates the second law of the thermodynamics in the metabolic networks helps to find the cyclic reactions (that should be eliminated) which led to infeasible fluxes (Pirt, 1965). Formulization of the thermodynamic constraints in the Thermodynamic based Metabolic Flux Analysis (TMFA) requires the Gibbs free energy values of reactions which could be determined by an estimation using group theory or from experimental data (Henry et al., 2007).
Thermodynamic constraints for determining the direction of reactions was firstly used in the reconstruction of the E. coli metabolic network. It should be noted that the thermodynamic constraints illustrate the feasibility and direction of the reactions (Kummel et al., 2006a, b; Henry et al., 2006). The reference concentration, 1 M, is considered for all internal metabolites in which the standard Gibbs free energy values of the reactions are calculated in this concentration. Since the real inside cellular concentration is about 1 mM, therefore, this concentration (1 M) does not demonstrate the real cellular concentration. The real values of the intracellular metabolite concentrations are between 5 μM≤C≤2 mM and also, it has to be noted that there is no any metabolite concentration from biochemical point of view if the intracellular metabolite concentrations are not at the mentioned range (Henry et al., 2007).
Interestingly, our method makes a correlation between the flux directions and metabolite concentrations. In this study, the effect of the metabolite concentration variations in the mentioned range on the Gibbs free energy values of the metabolic reactions and subsequently on the objective fluxes has investigated. The method has used in Genome-Scale metabolic network of Escherichia coli iAF1260. E coli iAF1260 metabolic network was able to more accurately predict the phenotypes of E coli MG1655 rather than the previous version of E. coli metabolic network reconstruction which has not involved the thermodynamic considerations. The biomass and the succinate production fluxes which have selected as the objective functions are simultaneously optimized using multi objective Genetic Algorithm in MATLAB environment. The calculations showed that the concentration and also Gibbs free energy changes could affect the objectives fluxes.
MATERIALS AND METHODS
Method Formulation
The study was carried out at the department of Chemical Engineering-Biotechnology,
Islamic Azad University, Science and Research Branch during 2007-2009. The reconstructed
networks could be resulted the BIGG (Biochemical, genomic and genetically) formation
data banks (Covert et al., 2003; Price
et al., 2003; Reed and Palsson, 2003). The
better understanding of the phenotype and also the objective functions in the
cell for development of the In-silico models which apply in the metabolic engineering,
is the goal of the computational models. As previously mentioned, the Flux Balance
Analysis (FBA) is a very good approximation for the estimation of the intracellular
behaviors. The FBA which could be obtained from the steady state mass balance
equations, is expressed as the following equation:
| (1) |
where, S is the stoichiometry matrix and V is the reactions fluxes. The numbers of reactions considerably are more than the number of metabolites, therefore, there are many possible solutions to solve the mass balance equations. In order to find a definite solution for problem, the constraint based analysis method used to confine the solution space. After imposing constraints, the metabolic network is optimized with respect to a certain objective functions. This optimization problem is a classical Linear Programming (LP) problem that could be solved using the simplex algorithm. Also, the system is limited by –α≤Vi≤β inequality, in which α and β are the lower and upper boundaries of the flux of each reaction. The constraints for the reversible and irreversible reactions are –α≤Vi≤β and 0≤Vi≤β, respectively. The constraint of the exchange reactions which lets the metabolite be discharged out of the cellular space is 0 = Vi≤β. According to the second law of the thermodynamics, the fluxes should move from the reactions with higher to the lower chemical potentials, which involved increased entropy (Maskow and von Stockar, 2005). Based on this law, many fluxes resulted from FBA are impossible. We have the chemical potential difference (Δμi) for each internal network flux (Vi). The second law of thermodynamic for each of these fluxes are defined as ViΔμi<0. As above mentioned, according to the second law of thermodynamics, the relation between the net direction of the biochemical reactions and Gibes free energy changes could be defined as follows:
| (2) |
where, Sgn() is the sign function, ΔGr is the Gibbs free energy changes and Vj is the flux of the reaction. The Gibbs free energy is obtained from the follow equation:
| (3) |
where, M is the active concentration (activity) of metabolite M, S and P are the substrate and the products, respectively, R is the global constant of Gas and T is the definite temperature of the reaction (298°K)
Also, there is a direct relation between the Gibbs standard free energy changes and the equilibrium constant of reaction as follow:
| (4) |
where, Keq is equilibrium constant of the reaction.
We reorganize the Eq. 3 as follows:
| (5) |
where,
The optimization procedure is expressed as follow:
Maximize Objectives
Subject to S'.V = 0
where, V is the vector of the fluxes, dj is binary parameter with the amounts of 0 and 1, α determines the upper boundaries and also S’ refers to the stoichiometry matrix, in which the columns of exchange reactions are eliminated.
By considering two above mentioned boundaries, it can be demonstrated easily that the Sgn(Vj) = -Sgn(ΔG) is always dominant (from second law of thermodynamic).
Model
Escherichia coli, the negative gram bacteria, is a significant choice
for reconstruction of the metabolic networks which is too controversial for
the metabolic studies and has the best specifications of a microorganism in
the field of genome, the objective function specifications and growth behavior
information. Since 1990 the reconstruction of the E. coil metabolic network
has been started and also is presently under investigation (Almaas
et al., 2004; Edwards and Palsson, 2000;
Flores et al., 1996; Covert
and Palsson, 2002; Price et al., 2004; Feist
et al., 2007). In the current study, the reconstructed genome scale
metabolic network of the E. coli iAF1260. This metabolic model also includes
thermodynamic considerations. The changes of Gibbs Free Formation Energy (GFFE)
and the Standard Gibbs Free Energy (SGFE) of reaction in E. coli iAF1260
has been already estimated for most of the metabolites and reactions (e.g.,
96% of metabolits and 84% reations). The GFFE and the SGFE values of the reactions
are calculated by the group contribution theory (Feist et
al., 2007).
The E. coli iAF1260 includes of 1972 metabolites and 2382 reactions. However, 300 exchange reactions are eliminated from the model and ultimately, the network includes of 2082 reactions and 1972 metabolites. Among the 2082 internal reactions of the network, there are near 81 lumped reactions in which the equilibrium state hypothesis has applied and also for which SGFE = 0. Since the concentrations of the metabolites can affect the SGFE and ultimately the range of feasibility of the reaction flux vector, the effects of metabolite concentration changes on E. coli iAF1260 metabolic network and the objective functions have investigated. The real values of the intracellular metabolite concentrations are between 5 μM≤C≤2 mM and also, it has to be noted that there is no any metabolite concentration from biochemical point of view if the intracellular metabolite concentrations are not at the mentioned range (Henry et al., 2007). In this work, by choosing a definite step size in this range, the effects of concentration on the fluxes and the objective functions are studied. The production of succinate has studied as a by-product in the TCA cycle of the E. coli (the wild type and the mutant types of this bacterium) are as objective functions alongside with maximizing of producing of the Biomass.
| Table 1: | Concentrations of nutrients in media |
Finally, the effects of the concentration changes in these two objective functions and in both wild type and the mutant types have showed. In this model, growth media contains Glucose, ammonium, sulfate, oxygen and phosphate have selected for the microorganism according to the biomass formation formula. The concentrations of these nutrients which obtained from experimental data are given in Table 1 (Henry et al., 2007).
THE optimization problem was solved using multi objective Genetic Algorithm accessed by the MATLAB (The MathWorks Inc., http://www.mathworks.com) modeling environment.
RESULTS AND DISCUSSION
The SGFE of the reaction of E. coli iAF1260 metabolic network have obtained from the data bases. For lumped reactions in which the Gibbs free energy cannot be calculated through the ordinary methods, it was assumed that the reactions are in the equilibrium state and their SGFE has set to zero (Henry et al., 2007; Feist et al., 2007; Mavrovouniotis, 1990, 1991; Maskow and von Stockar, 2005). The physiological concentration range for the intercellular metabolites has considered 5 μM to 2 mM. By selecting a specific step size 0.0002, the effects of the metabolite concentration on the objective functions fluxes of the wild type and the mutant types of E. coli iAF1260 (in which some reactions (genes) have been knocked out) have studied (Datta, 1992; Emmerling et al., 2002; Stols and Donnelly, 1997; Varma and Palsson, 1994). Selecting the biomass and succinate fluxes simultaneously as the objective functions and optimizing the model using the multi objective Genetic Algorithm, have applied. In order to study the changes of succinate yields in this range, five Mutant E. coli iAF1260 have selected and the reactions (genes) of the byproducts which have the least effect on the system biomass have knocked out. Succinate and biomass production of wild type E. coli and the mutants are listed in Table 2 in which the biomass and succinate production yields are provided on a basis of 10 mmol h-1 glucose and 1 g DW of cell. Our results represent in Fig. 1a-l which shows how biomass and succinate production fluxes of wild type and mutants of E. coli in Table 2 are changed versus concentration changes in 5 μM≤C≤2 mM range.
Figure 1 represents the effect of cellular concentrations changes on optimal growth as well as succinate production rate. As it is shown, until reaching a threshold concentration, the increase of cellular concentration yields the increase of the rate of both biomass and succinate production. It should be noted that by reaching to the threshold concentration, the increase of overall concentration of metabolites in the cell has no more effect on biomass and succinate production rates. This observation indicates that there is a definite cellular concentration in E. coli metabolic network, in which some cellular phenotypes like growth and succinate production rates are not sensitive to the cellular concentrations of metabolites.
| Fig. 1: | (a-l) Biomass and succinate fluxes changes for wild type and mutants types of E. coli versus concentration in 5 μM≤C≤2 mM range |
The obtained results demonstrated that the metabolites are in the average concentration of 1 mM in the cell. However, the results of this study indicated that by reaching to the threshold concentration, the metabolic reactions become insensitive to the cellular concentrations from thermodynamic point of view.
In addition, the obtained results revealed that there is a difference between the phenotype behaviors of wild and mutant types. When a gene (genes) is (are) knocked out from the genome of the organism like E. coli, the system shows different behavior before reaching to the threshold concentration. Interestingly, the obtained result in the current study confirms the previous studies which have been performed on the metabolic network of E. coli (as experimental work) by Palsson and co-workers which provided some mutants from E. coli, however, they have found that all mutants return into the biomass increasing phenotype after some generations as the cell metabolism objective function. They concluded that this finding demonstrates the stability of gene regulatory and transcriptomic network of E. coli in proportion to the imposed perturbation.
| Table 2: | Wild type E. coli iAF1260 and the selected mutants for succinate production |
| The reactions and corresponding enzymes for each knock outs are listed and the biomass and succinate production yields are provided on a basis of 10 mmol h-1 glucose and 1 gDW of cell | |
In addition, the obtained results showed that the mutant types of an organism not only tend to turn back into the wild type phenotype because of its transcriptomic stability, however, another factor is accessing to the thermodynamical stability which has evolved during the ages in the wild type. One of the most important factors in the evolutionary process of the mutant types is achieving thermodynamical stability in the non equilibrium conditions of living cells.
Also, the results showed that there are unsteady metabolite concentrations, in which we could analyze them from the metabolic engineering point of view. One of the problems of metabolic engineering is achieving types of the engineered strains which have fewer tendencies to return into the wild type behaviors. Obviously, if the organism returns into the wild type behavior, the efforts of the metabolic engineering will be wasted and the cell will use the sources in the ways other than the main aim of the strain. It is possible that the genetically engineered strains could not set their metabolite concentrations in the level of the steady state concentration. The findings resulted from the above diagrams demonstrate that the most serious result of this inability may be the problem of the cell in establishing the thermodynamical equilibrium state.
CONCLUSION
The major finding in the present study showed that there is a threshold cellular concentration in which before reaching to this concentration phenotypic behavior of cellular systems are not stable from thermodynamic point of view. This finding is a further emphasize on the significant of thermodynamic analysis of living systems in order to better understand the mechanisms of their functioning from the system engineering point of view. In other word in order to understand underlying interactions in the living systems and also how they affect the cellular behaviors, the study of biochemical networks and development of the experimental procedures and computational tools is so important.