From flux analysis to self contained cellular models

Abstract

Mathematical models for cellular systems have become more and more important for understanding the complex interplay between metabolism, signalling, and gene expression.In this manuscript, starting from the well-known flux balance analysis, tools and methods are summarised and illustrated by various examples that describe the way to models with kinetics for individual reactions steps that are finally self-contained. While flux analysis requires known (measured) input fluxes, self-contained (or self-sustained) models only get information on concentrations of environmental components. Kinetic reaction laws, feedback structures, and protein allocation then determine the temporal output of all intracellular metabolites and macromolecules. Emphasis is placed on (i) mass conservation, a crucial system property frequently overlooked in models incorporating cellular structures like macromolecular structures like proteins, RNA, and DNA, and (ii) thermodynamic constraints which further limit the solution space. Matlab Live Scripts are provided for all simulation studies shown and additional reading material is given in the appendix.

1 Introduction

Understanding complex systems, both in technical and non-technical sciences, relies heavily on mathematical modeling. This holds true for life sciences, particularly in systems biology, where mathematical modeling serves as a formalization tool for comprehending the intricacies of a system. Systems theory provides a framework for constructing models with various dimensions. One dimension pertains to the level of detail in the model, spanning from simple qualitative interaction networks to extensive, mass-conservative quantitative models detailing processes within cells and their fluctuating environments. Another dimension involves whether the model represents an average cell or individual cells within an environment. Modeling individual cells demands sophisticated approaches, such as employing population balance equations developed by Ramkrishna and colleagues (Ramkrishna, 2000) or adopting an ensemble modeling approach. Both methods necessitate distinct numerical schemes for resolution. Additional dimensions encompass whether the system is static or dynamic, and if the model requires structural elements, like an objective function, to explore the potentially infinite solution space.

The text at hand endeavors to initiate the analysis of a biochemical reaction network using flux balance analysis (FBA), a well-established method in cellular systems modeling. Therefore, only a basic understanding of reaction engineering is needed, while in-depth knowledge of microbial physiology commonly applied in biotechnology is unnecessary, since the examples are based on a simplified (toy) network. It concludes by constructing self-contained models that consider cellular macromolecular units linked to metabolism. Emphasis is placed on mass conservation, a crucial system property frequently overlooked in models incorporating cellular structures, leading to an incorrect degree of freedom for selecting system unknowns. While standard FBA typically considers a single biomass reaction, the proposed framework accommodates the structured nature of cells, highlighting resource allocation importance. Simple, self-contained models derived from this framework serve as a foundation for more intricate models.

A systematic procedure is presented, starting with a basic biochemical network devoid of macromolecular units. It delves into thermodynamic properties and introduces methods for a proper description avoiding inconsistencies with physical basics. Conforming to the dimensions mentioned earlier, the manuscript confines itself to detailed quantitative models for an average cell, analyzing static properties.

Therefore, it is not only suitable for beginners with basic background in mathematics (algebraic knowledge and differential calculus is required) and modeling but also for advanced students who could deepen their knowledge, especially when the macromolecular unit structure of the cell is introduced (here two different approaches are used).

After studying the text, the reader gains an understanding of the fundamental mathematical framework underlying models of cellular systems, as introduced in Section 2.1. Although a simplified (toy) model is used for illustration, comprehensive stoichiometric models of real cellular systems are available in specialized databases and can serve as a foundation for further exploration. Section 2.2 addresses basic thermodynamic challenges, enabling the reader to formulate systems of equations that avoid thermodynamic inconsistencies. Section 2.3 expands the discussion to include cellular models that incorporate macromolecular structures. A key outcome here is the development of objective functions that account not only for intracellular fluxes but also for the enzyme requirements necessary to sustain those fluxes. Section 2.4 introduces an extended equation system that integrates a larger set of macromolecular components and explores their potential roles in feedback mechanisms—such as how proteome allocation can influence intracellular flux distributions. Finally, the text introduces the reader to kinetic modeling approaches (Section 2.5). Whereas fluxes were previously treated as unknown variables, this section presents kinetic equations for individual enzymatic reactions. This allows for the analysis of the dynamic relationship between metabolic fluxes and intracellular metabolite concentrations, illustrated through the ongoing example of a small metabolic network.

Last but not least, hopefully lecturers and also well-advanced scientist will find interesting examples for their courses. A quantitative description always aims for a better understanding of the system at hand. Although the examples used are small, the presented methodology is general and can be applied to larger stoichiometric networks to determine intracellular flux distributions given measured uptake and production rates. Moreover, problems of resource allocation can be considered and will allow researchers to compare simulation results with own data. Larger differences should lead to the formulation and testing of new hypotheses that again will lead to new sets of experimental experiments. Additional reading material (books and standard publications) are provided in the Supplementary Appendix.

The examples provided adhere strictly to quantitative principles, employing physiologically meaningful parameters and standard units such as grams, moles, grams of dry weight , and hours for time constants. All simulations are executed using Matlab; live-scripts are provided for replicating the examples. Although software packages are available that support the reader by implementing models and performing simulations studies, only standard Matlab code is provided here. The advantage from the learning perspective is that the scripts can be directly used in the lecture or exercises. In typical courses, the focus is on the modeling approaches and how they are used to solve problems in systems biology and metabolic engineering. More details on available software tools can be found in the Supplementary Appendix.

2 Case studies

2.1 Mass balance in biochemical reaction systems

A common feature of all deterministic models that were developed for applications in systems biology, medical systems biology, or metabolic engineering is mass conservation. This is valid on the level of single bio-chemical reactions but also for larger networks to whole-cell models. For the studies presented here, mass conservation is obtained for single biochemical reactions which also result in mass conservation on the level of biomass synthesis; closely related to this, it allows a clear calculation of the specific growth rate if cell division is considered. The specific growth rate is defined as the change of biomass over time divided by the biomass (often in literature, a definition based on the concentration of biomass , that is, the ratio of biomass to the volume in the bioreactor system , is found that will, however, lead to inconsistencies when considering processes with a changing volume). With biomass, the weight of all cells in the experimental set-up is meant; the models developed in the text are based on a defined structure with metabolites and macromolecular units. To be mass conservative, the mass fraction of all metabolites and units (given in ) must add up to .

A standard notation with stoichiometric coefficients and reaction rate for a single reaction reads:

The numbers given in front of the metabolites are positive numbers (therefore they are written with value sign ) and describe the number of molecules involved in the reaction on a molar basis. They report on the left side the number of molecules that serves as substrates, and on the right side the number of molecules that serves as products. In the current form, these terms cannot be used for further calculations and must be translated into equations. As a general agreement, the stoichiometric coefficients on the left side are taken as negative values ( and ), while the coefficients on the right side are positive ( and ). The corresponding equation reads as follows:with the vector is the stoichiometric vector of the reaction . However, since the equation is only based on the numbers of the interacting molecules, we cannot infer, if the mass balance is closed. Here, the molecular weight of compounds come into play. A closed balance appears if we multiply the stoichiometric vector with the vector of the molecular weights of the metabolites:and as a results, a value of zero appears on the right side (scalar value since a row vector is multiplied with a column vector).

Since in a cell, a large number of reactions run simultaneously, the stoichiometric coefficients are stored in a matrix , called stoichiometric matrix (each reaction is represented in a column):

For compounds and reactions, the dimension of is ¹. Example 1

In this example, the notation for a chemical reaction is applied for a single reaction. A reaction equation reads (, , ):and with the molecular weights , , and , we see that Equation 3 is exactly fulfilled.

Variable is used to describe the velocity of the reaction. Important parameters that are introduced later on to describe this term are enzyme specific parameters like the turnover number , substrate and product binding constants , , and the equilibrium constant . Furthermore, the reaction rate also depends on the concentration of the reaction partners. In this way, could depend on a long list of concentrations and kinetic paramters .

Mass conservation on the level of single reactions implies that also for each compound in the system, an equation can be written down that sums up all processes that either increase or decrease the mass (or number of molecules) of the compound. Unfortunately, we cannot write down directly an equation that describes the course of the mass of the compound over time ; rather we can sum up all mass flow over time that changes the mass either in a positive or a negative way. The change of a compound over time is expressed with and the notation used in the engineering community is used in the following. This results in a ordinary differential equation (o.d.e.) for the mass of the compound or if we consider the molar number of compound . Typically, only biochemical conversion due to reactions is considered that includes uptake from the environment, excretion into the environment, or conversion into another compounds in the environment or inside the cell. Although diffusion is important in larger cells, it is not considered here.

The change of the mol number – that is proportional to its mass – of a compound is given by:

That is, the reaction velocity is weighted with stoichiometric coefficients and is multiplied with the total biomass . Note, that the sum is taken of all reactions with index for component with index (therefore, the stoichiometric coefficients always have two indexes). This is necessary, since, typically, a reaction velocity is given as a specific rate, that is, based on the total biomass with unit . In principle, we are done, and we could continue with the analysis of systems with a given number of reactions and compounds. However, in the current form on the left hand side, the change of numbers per time unit is given while on the right side reaction velocities appears that might depend on the concentration of the compounds . Therefore, a further step is needed here: To convert the left hand side, intracellular concentrations are defined by with unit , that is, based on the total biomass. Starting from and using the product rule for the derivative for , the final equation reads:with the specific growth rate that appears as dilution term which can be understood as follows: Assume that all rates that synthesize or degrade the compound are zero, the concentration of the compound will not stay constant but will go to zero due to cell growth (in every cell devision, the biomass is doubled). Equation 7 can be written for all compounds in the cell taking into account all reactions stored in the stoichiometric matrix . In this way, all stoichiometric coefficients for a single metabolite appear in the respective row of matrix .

Based on this fundamental equation, model analysis can be started. A standard tool here is the determination of all fluxes in a given metabolic network. Flux analysis focuses on balanced states, that is, it is assumed that all processes result in a quasi steady-state of the involved compounds. Furthermore, taking into account that in the upper equation the last summand has numeric values of different orders of magnitude, the standard equation for flux analysis simplifies:

Please note, that the dependencies of on the concentration vector is omitted here. In this way, the reaction rates are the unknowns of the equation system. Since the number of reactions rates (the column of the matrix ) is usually higher than the number of metabolites (the rows of the matrix ), the set of possible solutions of Equation 9 is infinite, and, usually is expressed in form of the null space of (null space vectors are orthogonal to each row, and therefore fulfill Equation 9 (see Supplementary Appendix). Since all linear combinations of null space vectors also represent solutions, some solutions are characterized by the following property: it is not possible to omit one of the reactions with a nonzero flux while for the remaining reactions, a different solution of Equation 9 can be found, or in other words, the solution cannot be decomposed. These special solutions are called elementary flux modes.

Example 2

For the bacterium Escherichia coli, the uptake rate for glucose in a minimal medium is estimated to be for a growth rate of . High values for intracellular concentrations for metabolites in central pathways are reported to be ², and . The dilution term is a factor 400 smaller than the uptake rate and it is justified to neglect it in this case.

Equation 9 will be the starting point for the analysis of our first network. One can assume that not all of the rates are unknown; typically good measurable rates are the carbon source, oxygen, , and nitrogen uptake. This allows us to split the rate vector into unknown rates (index ) and known rates (index ):and the number of unknowns is reduced. In general the number of rows corresponds to the number of compounds and the dimension of matrix is with reactions.

Formally, we are looking for a vector or that solves the Equation 10; this solution is called a flux distribution or a flux map. The number of solutions strongly depends on structural properties of matrix : in the cases considered here (for cellular networks, is valid), there are a infinite number of solutions³. Therefore, to find a physiological meaningful solution, one solution is selected that fulfils additional properties. One possibility is to choose a solution that maximizes or minimizes an objective function. Objective functions are chosen in such a way that physiological criteria are met, for example, bacterial systems tend to generate ATP as much as possible, grow as fast as possible or behave in a economic positive way in conditions with substrate surplus or scarcity (Schuetz et al., 2007)⁴. Mathematically, in this way, an optimization program is defined:with entries that select a single rate or combinations of rates to be extremal, and two “subject to” equations: constraints from the mass balance and boundaries for the fluxes.

2.2 Flux models without drain into biomass

On the first level, a simple network is considered, and the basic procedure to determine the flux maps is introduced. On a later stage, biomass synthesis will be considered. For the following examples, the stoichiometric matrix reads (the external components are not considered, therefore, there are only four rows and seven columns):and is illustrated in Figure 1.

FIGURE 1

Rates are numbered from (input flux) to . In this network, the number of reactions exceeds the number of components, and one can see that the rank (see Supplementary Appendix) of is and that the dimension of the null space (or kernel) of the matrix is . For the example, already a number of options are given, to select one flux distribution as optimal. In all cases, the input flux is fixed .

The calculation of flux distributions is shown for a number of different settings in Example 3. First, two different objective functions are applied and the different flux maps are compared Example 3a,b. For a selected objective function in Example 3c,it turns out, that some fluxes reach its maximal value. A closer inspection of the flux map reveals that thermodynamic laws are violated. After introducing the basic equations for a thermodynamic analysis for a single reaction and for biochemical reaction network, Example 3c is considered again and is now solved taking into account the presented equations.

Example 3

3a. Rate is selected to be maximized. 3b. Sum of all fluxes is minimized. 3c. Rate is maximized. Optimal solutions for seven different settings are shown in Figure 2. The flux maps differ significantly; in case 3a, all of the incoming flux into compound could be converted into . In case 3b, we note, that minimizing the sum of all fluxes, results in only two fluxes that were active while all other fluxes are zero. From all possible elementary flux modes for this system (shown in the Supplementary Appendix), the calculated flux distribution has overall the minimal sum.

FIGURE 2

For 3c we note that–independent from the input flux–the objective always shows the same value. This value corresponds to the upper limit that was set by the boundaries. Boundaries limit the solution space and are based on (experts) knowledge on the system; they could be positive or negative, depending on the problem formulation. Since the direction of the fluxes are not known before in general, positive and negative values are used. How can we interpret this result? A closer and detailed look in this case reveals that not only flux shows a special behaviour but also fluxes and . Both change their values according to match Equation 10. This means that the fluxes are running in a cycle, represented by the three rates shown in Figure 1.

2.2.1 Thermodynamics

Example 4

For a reaction with one substrate and two products (cf. Example 1):

The forward and the backward reaction are given by , , if simple kinetic rate laws (mass action) are applied. Here the reaction rate is proportional to the involved metabolites. Setting forward rate and backward rate equal, one get:

In this case, Equation 13 reads:and one obtains finally by combining the last two equations:

If , that is, , the standard value is negative (the log is smaller than zero for values smaller than 1), indicating that under standard conditions, the reaction would take place from left to right. Reaction kinetics could be more sophisticated in case that enzymes act as catalysts (see below).

Taking advantage of the logarithmic function, the exponents can be written before the symbol; furthermore, using matrix operations, Equation 13 can be re-written in a scaled nomenclature for all internal reactions of the network using the stoichiometric matrix and vector representing the logarithmic values of the concentrations :with is the vector with entries of the scaled free Gibbs energies of the reactions (that is divided by the term ). Only internal reactions are taken into account, since information on external concentrations of metabolites are not considered here.

From what we have seen before, a simple relationship between fluxs and values is obtained:

That is, the product of the free Gibbs energy and the respective flux is always negative. For more than one reaction, the following notation is used:

2.2.2 Thermodynamics in networks

At this point, we have to note, that the values of and therefore also the equilibrium constants are not independent in a network. If we consider again the cycle with reactions , , we see two ways from to ; the first one directly via and the second via and (Figure 1). In other words, one can think on water that falls down a fall on two ways, however, at the bottom the energy status for both ways is equal and it is not possible to bring back the water on top without bringing in energy.

In the last section, we have seen that fluxes are coupled with metabolic concentrations via the Gibbs free energy. Considering large networks, this would require additional state variables (either the Gibbs free energy for reaction , the chemical potential for component , or the concentration itself) to obtain a valid flux distribution. In the following different ideas are introduced to guarantee that the fluxes are not running in a cycle. The meaning of cycle is strictly limited here to thermodynamic considerations; other types of cycles such as futile cycles or control loops which are thermodynamically feasible are not meant.

As a starting point, an equilibrium is considered, that is, the net rate of all fluxes is zero as well as the three values. From above, the following set of equations is obtained (the equations are valid only in equilibrium, here we omit index ):and one can see by close inspection:

The equation can be re-written with the equilibrium constants of the reactions (Equation 18):

Can we obtain these equations also in a more systematic way? The cycle given by the three reactions can be obtained by the following equation; if only internal reactions are considered, that is, there is no exchange of mass with the environment, then, the null space of matrix represents the reactions of the cycles in the network. In our case, the matrix of internal reaction reads:and one obtains for the null space:

Finally, we see thatis exactly Equation 23.

The last equation is not only true for the equilibrium. The Gibbs free energy can be written as a function of the stoichiometric matrix and the vector of the chemical potential :

Starting from the observation from above that the null space vector of contains all cycles in the network, the following relationship is valid:

A multiplication of the last equation from the left side with the vector gives:

That is, the values for must fulfill additional constraints; a multiplication with all vectors in the null space of must be zero.

Example 5

Considering the loop in the network above, one gets:

Now, assume that in the solution flux is zero, that is ; as a result, we get the constraint from the last equation:

In this case, a solution with a flux from to via reaction and with positive values is forbidden; if there is no flux between and due to a reaction equilibrium, this is also true for every different way from to .

To incorporate these conditions in the linear framework used so far, we take advantage of a method described by Schellenberger et al. (2011). The methods uses integer variables for each reaction/Gibbs energy with values only or . One difficulty given by the equations above is, that a vector with zero entries for the values is always a solution to Equation 30 which is, however, undesired. To overcome this, the following equation system excludes value equal zero for the vector:with is a small number while is a sufficiently large number. For the two cases (i) we obtain:

That is, for negative flux rates, the corresponding Gibbs energy is positive. In case (ii) we obtain:

That is, for positive flux rates, the corresponding Gibbs energy is negative. This set of equation can now be written in matrix form in a very compact way. The full vector of variables that must be solved, now is as follows: and the equation system with equality and inequality condition reads (sub-matrices , and 0, as well as , and have appropriate dimensions):

Note that matrix is used to make clear that only internal reactions of the network are considered for the thermodynamic analysis, that is, all external rates have entry zero in . Matrix is the null space of the internal reactions as used above. Note, that cases with zero fluxes in the cycles might lead to inconsistencies because the corresponding value cannot be zero. However, in the Supplementary Appendix SA method is described that can be used for a verification, if a given flux map with zero entries in the internal flux vector is feasible.

A drawback of the procedure is, that additional variables including integer variables are needed that generate challenges while solving the equation system. An elegant way, described in (Beard et al., 2004), to avoid that fluxes are running in cycles is to formulate a non-linear equation based only on the sign of the fluxes (vector ) and the sign of the cycles (vector ) for :

Simply speaking, if a flux solution runs in a cycle, the left and the right hand side are equal; if this is not the case, then the left hand side is always smaller than the right hand side. Please note, however, that the strong constraint given in Equation 30 are not valid here; if there is potential between two components, a flux might be zero. In our case, for the sign vector we get . In case that the flux distribution follows not the cycle, for example, , we get . In the normal case that the null space is high dimensional, the given equation should be applied to all cycles in the network. However, to determine the complete set of cycles is difficult, since linear combinations of null space vectors are also in the null space.

A further alternative, including the concentration of the metabolites and replacing the Gibbs energy by equations give above, allows us to find a further set of equations that could be used: Above, the condition (Equation 21) was given that relates the sign of the Gibbs energy to the sign of the rate. Now this is applied to the cycle in the network. We restrict the analysis to a given flux vector, however, we only use the sign of the vector; . The set of equation read:

Example 6

The set of equation is applied for the cycle in the network, and let us assume that and are positive and is negative, that is, fluxes are running in the cycle. Equation system (Equation 22) can be copied with slight modifications (we directly use the equilibrium constants and replace the equal signs with and , respectively):

From the first two lines, it follows:with the condition for the equilibrium constants , the last line from above reads:

That is a contradiction to the first two lines; the flux cannot run in the cycle. The presented equations are sufficient for a loopless network.

Example 3c (re-opened)

For the example, an extended optimization problem is formulated taking the concentrations of the compounds into consideration. The general form now reads:with takes into account only internal reaction rates and no exchange reactions. As above, algorithms used in the study require upper and lower bounds that are given in the last two lines. Example 3c is solved with equation system (Equation 42), and we get a different flux map fulfilling all constraints now (Figure 3).

FIGURE 3

However, there is a broad range for the values for the concentrations fulfilling the given constraints. In analogy of selection one extraordinary flux for the optimization procedure, we proceed the same way and select from all possible values for a set of concentrations one which fulfills additional specifications. As proposed in literature, for a given flux map, all (absolute) values for the internal reactions are determined (since the direction of a reaction is open, we are interested only in the quantity of ), and one tries to maximize the minimal value of these (Max-min driving force (MDF) design principle (Noor et al., 2014)). Formally, a different optimization problem is formulated based on given flux map for the internal reactions only:

Example 3c (continued)

Figure 3 left shows the updated and thermodynamically valid flux distribution. In the middle, shown in red, are the values for the concentration of the compounds without and with the MDF design principle. Right: The minimal value for the negative is higher in the MDF case. An interesting observation is that, although is higher than the other values, the corresponding flux is nearly zero. We note that, so far, fluxes and concentrations are only weekly coupled. Later on, this problem can be resolved by taking into account that fluxes are proportional to enzyme concentrations, and complete enzyme kinetics can be exploited.

To complete the example, we perform a multi objective approach to maximize flux , and minimize the sum of fluxes. The results of such a calculation can be plotted as a Pareto front. The values on the front are characterized by the following property: when we move along the front, only one of the two criteria is improved while the second one is impaired. For point that are not on the front, we can improve both criteria. The arrows in Figure 4 left indicates the direction for the improvement in both directions; crossing the front is not possible. Also, the values for the four internal reactions show a linear behavior.

FIGURE 4

Lessons learned from Example 3: To determine valid flux distributions for cellular reaction networks, mass conservation as well as thermodynamic principles must be considered. Depending on the choice of the objective function, flux maps may differ. The selection of the objective function depend on the experimental set-up. For growth of bacteria in surplus carbohydrate the maximization of the growth rate is an appropriate choice. However, in case that substrates are limited or for applications in metabolic engineering, different objective functions are more suitable.

2.3 Flux models with drain into biomass

In this section, the network is extended to take into account that metabolites in the network often serve as precursors for biomass synthesis. Two different variants will be considered as shown in Figure 5; first, the standard notation used in FBA is introduced, then, a modeling approach with several macromolecular units such as protein, RNA and DNA is considered.

FIGURE 5

A standard approach in FBA adds a single reaction, also called a pseudo-reaction, representing 1 g of biomass. In this single reaction, the complete metabolism for biomass assembly is summarized. For the example network from above, we use on a gram basis:

The equation is converted with the molecular weight of the compounds and the stoichiometric matrix has to extended with an additional column. On a molar basis, the equation reads in the FBA notation:

In FBA, the right side in the last equation remains empty, that is, biomass is not explicitly taken into account as compound in the stoichiometric matrix and therefore, no additional row is necessary. This can also be seen in Figure 5 (right side); the reaction arrow points out of the system.

Aside note:

It is worth to have a closer look on the last equation since we changed the units from given typically in to the growth rate , given in . From our general approach for mass balance equations, we derived Equation 8; applying this to a compound (concentration is ) that is only synthesized in one reaction but not degraded (for examples macromolecules (with unit ) or a representative for the overall biomass), one obtains:

Normally, as in the examples given in the beginning, the stoichiometric coefficients are dimensionless (essentially they are ). For the special case of the biomass reaction, we write for example:but we do not consider compound biomass as an additional state variable. The o.d.e. for compound is used and is now transformed by excluding the biomass concentration from the expression of and writing it directly to the stoichiometric coefficient:

The last equation is re-written:

We note a change of units in the stoichiometric coefficients; is that is now changed for in . If, like in our case, the reaction equation describes the synthesis of the overall biomass, the rate of synthesis can simply be replaced by the growth rate , and the stoichiometric coefficients with unit represent with biomass scaled values.

In Example 7 the drain of metabolites to generate biomass is considered and as in the previous case, two different objective functions are compared. In addition, constraints by the available amount of enzyme for the network is studied, resulting in four different scenarios.

Example 7

The network considered so far is extended by a single biomass reaction (Equation 44). Simulations results are shown in Figure 6 with two cases: (i) maximization of growth rate, and (ii) minimization of the sum of all fluxes.

FIGURE 6

For the next level, a further extension of the stoichiometric/thermodynamic approach is introduced. It becomes obvious during the last decade, that the standard approach has a strong disadvantage; the biomass constituents are not considered, and therefore, any feedback from the biomass parts back to the metabolic network is missing. This is true especially for proteins like enzymes that are involved in any metabolic processes. Measurement of the proteome now offers several possibilities to extend the flux analysis approach by taking into account not only different protein classes or sectors, but also information on individual enzymes. In a step-by-step approach we provide, and also combine, several methods from literature. The intention is to show that protein allocation has a strong impact on flux distributions, which has implications for the understanding of cellular physiology and but also is important for applications in Metabolic Engineering where one is interested to redirect fluxes to desired products.

Considering individual enzymes, from kinetics, it is convenient to split the mathematical term into two parts in a first step. Such rate laws depend on the concentration of all involved substrates , products , and the overall amount of enzyme in the reaction assay. If the enzyme is under control by inhibition or activation, also the concentration of the effectors come into play. The kinetics can be written in the following form with a set of parameter :

From here we see, that because for the term , the following relation holds: for a irreversibel reaction. Moreover, it is observed that the total amount of protein during different conditions does not change too much, indicating that the sum over all enzymes could (i) be constraint, or could (ii) be subject of optimization, since a minimal sum of metabolic enzyme could allow the cell to spend more proteins in ribosomes to accelerate growth. Both ideas ((i) and (ii)) are described in literature as GECKO (Sanchez et al., 2017) (we used a slightly modified version here) or ECM [enzyme cost minimization (Noor et al., 2016)].

Taking into account additional variable for the single enzyme (for simplicity we omit the zero for the total amount of this enzyme), Equation system (Equation 11) is complemented. With the same objective function as before, one gets for case (i):with is a diagonal matrix with values , corresponding to the unknown rates ; is a parameter that limits the total amount of enzyme that is available for the network. Since the direction of the fluxes is open we must use the absolute values of the rates here.

In case (ii) a different objective function is used; namely, the sum of all enzymes in the network is minimized. Here, the equation system reads with weighting factors, for example the molecular weight:

In case that only positive fluxes are considered (this can be achieved by splitting reversible fluxes into two positive ones), the two inequalities can be re-written in matrix notation with matrices/vectors with apparent size:

Example 8

For the given network, a rough estimation of the can be done. From published proteome data, a mean value of transporter molecules per cell is obtained. For a calculation with a standard uptake rate in the range of , the number of molecules per cell must be transformed into the unit . Here, we need Avogadro’s number and an average cell mass :and one obtains for the :

Example 7 (continued)

Results of the different simulation studies introduced are shown in Figure 6. For a fair comparison, the growth rate was fixed in case of minimizing the sum of fluxes or minimizing the sum of enzymes to the same value that was obtained while maximizing the growth rate. The weights from above are taken as one. As can be seen, depending on the objective function, the flux distribution changes for nearly all rates ( and were set before and do also not change). The right plot in the Figure shows the sum of all fluxes. As expected the second objective function results in lower values, but the difference is only about 16% with respect to the first objective.

Lessons learned from Example 7: The choice of the objective function and additional constraints show a strong impact on the flux distribution. However, while some fluxes does not change (for example , , and ) other fluxes show drastic changes.

2.4 Flux models with several macromolecular units

The approaches so far still have a couple of limitations. The results are scalable since, typically, the main incoming fluxes must be given, and cannot get as a result from the optimization procedure; second, the biomass equation is an approximation, and does not account for different needs of the cell for growth and maintenance. In the next step, we omit the single biomass reaction and replace it with macromolecular units that are not subject to degradation (see Figure 5 right).

From the basic mass balance Equation 8, a steady state can only be reached, if the dilution term is considered. Also, to keep our mass conservative approach alive, the growth rate must be the result of incoming mass flow minus outgoing mass flow. Fortunately, from the basic equation, a simple, but powerful relationship could be derived for the specific growth rate (Doan et al., 2022):

The multiplication of from the left side with the molecular masses always results in entries 0 when the respective reaction is completely mass conservative. All entries that are not zero contribute either positively or negatively to the specific growth rate . The equation can directly plugged in our standard equation, and one obtains an extended steady state condition:

In comparison with standard FBA, an additional matrix comes into play that takes into account the concentration of all compounds and the macromolecular units. Fortunately, the system is still linear (if thermodynamics is not considered) in the fluxes and can be solved with standard tools. In addition to the flux distribution that were obtained so far, additional solutions appear that describe fluxes into the biomass. However, a close inspection of the equation reveals that the composition of the cell must be known (represented in vector ). Rough values for the macromolecular composition can be found in literature, but concentrations values for small molecules are difficult to obtain. However, for first calculations, we can set these values to zero.

In Example 9, a simple model structure is used to analyze the influence of the biomass macromolecular structure. It is assumed that these elements are not actively degraded, and therefore only diluted by biomass growth. This allows model reduction, and, as shown in the extended case, also allows to study the influence of feedback structures, that is, the macromolecules influence their own synthesis.

Example 9

A simple scheme with one uptake reaction, one reaction that excreeds the metabolite directly back into the medium (this is based on experimental observations and is termed “overflow”⁵), and two cellular units and is considered. Metabolite , a catabolic product, represents a compound from central metabolism or an amino acid that is needed to build macromolecular units like protein, lipids, or DNA.

The stoichiometric matrix for the example shown in Figure 7, left, reads as follows, taking into account that a higher number of monomers are needed for 1 mol of a macromolecular unit (that is, and have large values):with the first column represents the input flux as before. The molecular weights are in this case and . However, the quantities are not independent; for example, we get for the vector with the molecular weights:

FIGURE 7

While the masses of the compounds only fulfill the conservation condition (to sum up to biomass). A closer look at the o.d.e’s for the biomass components reveals a further interesting property that can be used for the numerical solution. For macromolecular units and , the respective equations read:

One gets for the steady state:

The nice thing with the last equations is, that scaling of the rate of synthesis with the molecular mass fractions always result in the specific growth rate, that is, the rates of synthesis are not independent (here: ). Vice versa, knowing the specific growth rate, that in the example can easily be determined:

The rates of synthesis for the macromolecular units can be obtained, and does not require the solving of an additional equation system. For this model, the null space can be used to determine all solutions of the system, and for the example, we obtain for the special case that the mass fraction of compound is negligible :

The first null space vector shows the way of the substrate through the network into the excreted overflow product (column 1) while the second vector represents the growth mode of the network (column 2), that is, the way from the substrate into the macromolecular units. A close inspection of the third element in the second vectors reveals that the element can be rewritten as , that confirms the statement from above.

As indicated in Figure 7a reduced model can be obtained with only one single flux into the structured biomass (2 fractions). In the reduced model, the biomass is structured in two parts with only one rate.

Aside note:

If the stoichiometric is re-written for the two biomass reactions in such a way that the mass composition is directly included, then the null space directly shows that all rates of biomass (or macromolecular) synthesis are equal to one for the growth mode and for the assumption . For the example above, then reads:and the product is:and the null space is:

Example 9 (continued)

Here, the reduced model is considered and a feedback structure is introduced as shown in Figure 8; both fractions of the biomass are needed to sustain the incoming flux (fraction ) and drain to biomass flux (fraction ). Here, we note, that the allocation of both fractions play an important role for the system. Depending on the incoming flux, the biomass must be distributed in such a way that both rates can sustain their fluxes. For a first simulation, the biomass composition is varied and the concentration of metabolite is set to zero. Since the total biomass is scaled to , only one control variable (here fraction ) is needed to explore the complete solution space shown in Figure 8 right.

FIGURE 8

The equation system starts from the standard equation system (Equation 58). Now, the rate vector (or parts of it) is replaced by a simple kinetic expression , that is in our case, , and (in this way, the system is still linear):

However, we must add the conservation equation, saying that the sum of and is constant and represents the overall biomass concentration :and we optimize for the specific growth rate .

Figure 8 shows the network structure and the simulation results. Since only only degree of freedom exists, the complete (linear) solution space could be determined (incoming flux , , and as well as the specific growth rate in dependence on fraction ). However, the results are not intuitive. The optimal value for the growth rate is shown with a red square. Smaller values than the optimal one for (values on the left side) result in a higher growth rate that could not be realized by the incoming flux (note that is negative). Higher values of (values on the right side) lead to a lower growth rate because is too high and fraction , responsible for its own generation, is too small. With decreasing rate , rate must increase further to reach a steady state.

In the next Figure 9, we show two simulation studies with this simple model. First, we mimic a limitation of substrate in the medium (two left plots). This is realized by changing parameter and the x-axis in this case represents the limiting substrate concentration (arbitrary units). The plots on the right side show a simulation study when we assume that not all is available for or . This case mimics for example, the synthesis of a heterologous protein.

FIGURE 9

The examples already demonstrate a self contained model system. In comparison to the previous models, substrate intake is limited by the available amount on sector . Actually, increasing the input flux is possible by shifting the resources to , however, of costs of synthesis of biomass. The optimal growth rate is a compromise between substrate uptake and biomass synthesis that, of course, depends on the kinetic and stoichiometric parameters used. The model also reproduces the behavior observed for protein sectors in Escherichia coli in the two applications. Studies on proteome data for different substrates and growth conditions reveal that the fraction for carbohydrate transport and central pathways is decreasing with the growth rate while the ribosomal fraction is increasing Schmidt et al. (2016). For this case, one also could expect a different behavior of the system with respect to resource allocation. The optimal solution indicates that the loss of resources due to the additional protein production is equally distributed to both fractions, however, a different solution could be that one sector stays constant while the load is given to only the second sector. One interesting observation is, that in contrast to literature, the behavior of the growth rate in the first application is not linear but shows a hyperbolic outcome and approaches a threshold. This is understandable since - although substrate uptake is linear and not saturated–the resources can be distributed either to the incoming uptake unit or to the synthesis unit.

The last set of examples introduces a new aspect in modeling cellular systems, namely, the dependencies of the reaction rate from the concentration of substrates, products and parameters. The analysis start with the network from Example 3 and then extends it with kinetics. The final case takes into account that the resources to generate the biomass macromolecular structure is limited by the proteins. In this way, a self contained cellular model is obtained.

Example 10a

In the next step, our example with 4 metabolites is considered again and extended with three macromolecular units. In Figure 10 the structure of model variant is shown that is used for three different cases.

FIGURE 10

In the first variant, the macromolecular units are fixed and the stoichiometric matrix contains the drain from the metabolites into the respective units. It is assumed that the sum of the mass fraction of the metabolites (, , and ) is that is 10% of the entire mass. This variant is yet not self containing; therefore the input flux must be fixed. Here, two objectives are considered: maximizing growth rate and minimizing the sum of fluxes (by a given input flux). Results are shown later on and are compared with the other cases (see Figure 12).

2.5 Kinetics–the step to dynamic models

Now we go back to Equation 51 and fill the kinetics with “life”, that is, we consider the dependencies of rates not only from enzymes, but also from the concentration of reaction partners. In literature, numerous kinetic expression can be found taking into account the number of catalytic centers as well as the number of substrates and/or products. The best known kinetics is the Michaelis-Menten (MM) kinetics that reads:with the substrate concentration and the half saturation parameter (also called Michaels-Menten constant). This kinetic describes the irreversible conversion of one molecule of substrate into one molecule . To be more flexible, for the following studies, a reversible reaction is used here, that allow a general conversion of molecules of into a number molecules of products (Wortel et al., 2014):with the abbreviations and in the denominator; and are the respective binding constants for substrate and product. The equation can be structured into three parts: one capacity term indicating the maximal value of the velocity, a saturation term between zero and one, taking into account the number and type of binding centers for substrate, product and possible effectors (not considered here), and finally, a thermodynamic term that describes the capability of reaction reversibility (Noor et al., 2013):

As seen above for a simple kinetic, a relationship in case of a reaction equilibrium could be obtained by setting . Here we obtain:and the equilibrium constant is obtained like before:

Please note, that in this framework, the value of determines the sign of the reaction (for example, if , the reaction rate is always positive).

Aside note:

In reversible enzyme kinetics, the equilibrium constant plays an important role and determines not only the direction of the reaction, but also the Gibbs energy. The reaction kinetics for a simple reversible reactions reads for :

We note that the kinetic parameters of substrate binding and product binding as well as the reactions parameter for the forward and the backward reaction and , all contribute to the overall equilibrium constant ; this relationship is called the Haldane relationship and can be written:

Keeping and constant and varying in the order of one magnitude, the Gibbs energy changes (left plot in Figure 11, and ). This can be seen by a simple calculation; from the expression for , the difference for two values of with a factor of 10 is taken:

FIGURE 11

However, taking Equation 75 to determine the amount of enzyme needed for varying product, a highly sensitive behaviour is observed (right plot in Figure 11). Near the equilibrium constant, the amount of enzyme needed to maintain a flux of the amount of enzymes raises very fast.

Furthermore, the above considerations on thermodynamic properties in cycles are taken into account, that is, the equilibrium constants are not independent, but must fulfill the condition given in Equation 24.

Example 10b (now with kinetics)

Now the same network as in Example 10a is considered, but kinetics for the metabolic network are taken into account. The system is now self-contained with a minimal number of constraints: mass fraction of metabolites is again fixed and the concentration of the unitss are given. For the models considered so far, the concentration of metabolites only had a modest role; they were only important for thermodynamic considerations. Figure 12 shows the results. Please note that for the equation for one macromolecular unit, the following equation holds in steady state:

FIGURE 12

Since in this case, the concentration of the unit is fixed and is only calculated via the metabolic network, the rate of synthesis can be directly obtained from the equation, and no further dependencies, for example from the metabolites, are necessary.

Example 10c (resource limitation)

In the last case, the first unit is responsible for all enzymes that are active in the metabolic network, that is contributes to (there is a fixed amount that is increased by the sum of enzymes). Unit is in charge for the synthesis of and for its own synthesis. Unit is is obtained by closing the mass balance, accordingly, its rate of synthesis can be determined from its concentration. Still, the sum of the metabolites corresponds to 10% of the entire cell mass (a minimum concentration is also required, that is, the concentration of a metabolite cannot be zero). Figure 12 compares results for three cases. For the equation for a self controlled unit (example 10c) a different equation holds in steady state:

A simple first order dependency in the form would immediately result in and arbitrary, which is independent from the complete network and from the environmental conditions. Since this is not useful, a simple kinetics, coupling the metabolic network to the synthesis unit, is used here:

The optimization problem is now as follows (for all cases in Example 10):

The substrate uptake rates in all cases are comparable (left plot) while the flux maps differ: rate is used in all but one case while rates are used in the positive direction only in the first case. For example 10b and c, enzyme is allocated for the uptake rate and for the synthesis of metabolite which fill up all other metabolites in the reversed direction (rates are negative).

The last model variant finally is used to mimic the behavior for a decreasing substrate concentration. In the previous example, the enzyme for substrate uptake is working in saturation, that is, the concentration of the substrate is far larger than the affinity constant and the flux can be approximated by a linear relationship . Now, the substrate concentration is changed and the results are plotted over the concentration of .

Although the input flux is linear, the resulting growth rate tends to saturate for high values. This is based on (i) the saturation kinetics of the individual enzymes in the network, and (ii) on the limited capacity of the enzyme fraction that can be allocated for substrate uptake; please note, that for the synthesis of enzymes (and unit ), unit must be synthesized too. Therefore, as can be seen, a typical behavior for the course of the units can be detected; unit which is self-controlled increases with the substrate (and with growth rate ) while unit (remainder) decreases. The amount of enzyme needed is included in but is not visible in the simulation. For high values of the substrate, only the enzyme for substrate uptake and contributes significantly; for smaller values, the pattern changes, and enzyme is needed to maintain the flux from to and which was not necessary before. Interestingly, the yield (that is biomass generated from substrate) increases for small values of ; this is due the lower flux via that besides generates also a by-product (see Figure 10).

We finally compare and illustrate the results of the simulation studies with biomass using elementary flux modes (see Supplementary Appendix for a definition and an example). For both networks, elementary flux modes with substrate uptake (from ) and building biomass are considered. Figure 13 shows the relationship between the obtained growth rate and the sum of fluxes (only metabolic fluxes). The black line indicates the maximal growth rate that can be obtained from substrate uptake of of .

FIGURE 13

As can be seen, the macromolecular unit model has more elementary flux modes than the single biomass model, showing a higher flexibility in providing precursors for biomass. However, the elementary flux mode with a higher growth rate than (the red symbol on the right side of the black line) includes an uptake of a second substrate (that is, the corresponding product is not build but serves as a substrate). As expected, while minimizing the sum of all fluxes, an elementary flux mode is directly obtained, the values obtained from example 10b and 10c (only growth rate maximization) are only near an elementary flux mode but do hit them perfectly.

A further analysis can be done by focusing on one node in the network and characterizing the incoming fluxes, that are named supply fluxes, and the outgoing fluxes, named demand fluxes. In steady-state, supply and demand are exactly balanced, that is, the sum of supply fluxes must match the demand fluxes. Moreover, since kinetic expressions are used for the last examples, the influence of metabolic concentrations can be studied. In the last example shown in Figure 14, the distribution of enzyme over the substrate input abruptly changes. In the following figure, two values of substrate input, corresponding to growth rates (high) and (low) are compared.

FIGURE 14

In the network, the demand is split in several parts. Metabolite is the center of the network because it is the first metabolite in the network. Hence, it must feed all other metabolites as well the direct drain to biomass. In the left plot in Figure 15, the demand bars show that main part of the incoming fluxes are used for the network demands (yellow and red part) while the flux directly into the biomass (magenta) is smaller. Since the concentration is small in comparison to the units, the dilution term (green) is very small. To understand, why in the low mode a shift of the allocation of enzymes from to is observed (Figure 14), a look on the concentrations of and reveals also a change here. Since reaction is reversible, in the high mode, the situation is not favorable enough because is lower (and is higher); in contrast, in the low mode, is higher (and is lower)that prevent that alone feed the other metabolites in the network and, consequently, comes into play to feed the network “from the other side” (see Figure 16, left plot). In this way, for example, the demand on metabolite is satisfied via two complete different ways; in the high mode with and the reverse reaction and in the low mode directly via . The dilemma is also illustrated with the enzyme allocation plot in Figure 16, right. For a low concentration outside, more enzyme must be spend for substrate uptake (compare bars for ). However, the complete amount of enzyme is restricted, and for reaction , enzyme is not available in high amounts. Consequently, enzyme spend to reaction is reduced and enzyme for reaction is significantly increased in the low mode. The example also demonstrates the high sensitivity of the system; while the change in the concentration of is marginal ( 7%) it goes along with a dramatic change in the proteome pattern.

FIGURE 15

FIGURE 16

3 Summary

The provided examples demonstrate an increasing complexity starting from a detached network with a couple of reactions and ending with a prototype of a cellular system with metabolism and macromolecular synthesis as well as kinetics for the metabolic network. While the simple network structure without biomass reactions must be scaled by the incoming flux, the last examples provide self-contained models. Thermodynamic considerations play an important role for the outcome of the simulations, but could be considered in the kinetic reaction expressions directly, bridging pure fluxes with metabolite concentrations.

Starting from flux balance analysis, using different objective functions to finally a three macromolecular unit model, increasing levels of complexity are developed and demonstrated by numerical examples. A final summary of basic equations and examples are given in Figure 17.

FIGURE 17

References

• 1

  BeardD. A.BabsonE.CurtisE.QianH. (2004). Thermodynamic constraints for biochemical networks. J. Theor. Biol.228 (3), 327–333. ISSN 0022-5193. 10.1016/j.jtbi.2004.01.008

  • CrossRef
  • Google Scholar

• 2

  DoanD. T.HoangM. D.HeinsA.-L.KremlingA. (2022). Applications of coarse-grained models in metabolic engineering. Front. Mol. Biosci.9, 806213. 10.3389/fmolb.2022.806213

  • CrossRef
  • Google Scholar

• 3

  EbrahimA.LermanJ. A.PalssonB. O.HydukeD. R. (2013). COBRApy: constraints-based reconstruction and analysis for python. BMC Syst. Biol.7, 74. 10.1186/1752-0509-7-74

  • CrossRef
  • Google Scholar

• 4

  HeirendtL.ArreckxS.PfauT.MendozaS. N.RichelleA.HeinkenA.et al (2019). Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat. Protoc.14, 639–702. 10.1038/s41596-018-0098-2

  • CrossRef
  • Google Scholar

• 5

  KremlingA. (2021). A counting-strategy together with a spatial structured model describes RNA polymerase and ribosome availability in Escherichia coli. Metab. Eng.67, 145–152. 10.1016/j.ymben.2021.06.006

  • CrossRef
  • Google Scholar

• 6

  NoorE.Bar-EvenA.FlamholzA.ReznikE.LiebermeisterW.MiloR. (2014). Pathway thermodynamics highlights kinetic obstacles in central metabolism. PLOS Comput. Biol.10 (2), e1003483–12. 10.1371/journal.pcbi.1003483

  • CrossRef
  • Google Scholar

• 7

  NoorE.FlamholzA.Bar-EvenA.DavidiD.MiloR.LiebermeisterW. (2016). The protein cost of metabolic fluxes: prediction from enzymatic rate laws and cost minimization. Plos Comp. Biol.12, e1005167. 10.1371/journal.pcbi.1005167

  • CrossRef
  • Google Scholar

• 8

  NoorE.FlamholzA.LiebermeisterW.Bar-EvenA.MiloR. (2013). A note on the kinetics of enzyme action: a decomposition that highlights thermodynamic effects. FEBS Lett.587 (17), 2772–2777. 10.1016/j.febslet.2013.07.028

  • CrossRef
  • Google Scholar

• 9

  RamkrishnaD. (2000). Population balances. Academic Press.

  • Google Scholar

• 10

  SanchezB.ZhangC.NilssonA.LahtveeP.-J.KerkhovenE. J.NielsenJ. (2017). Improving the phenotype predictions of a yeast genome-scale metabolic model by incorporating enzymatic constraints. Mol. Sys.Biol.13 (8), 935. 10.15252/msb.20167411

  • CrossRef
  • Google Scholar

• 11

  SchellenbergerJ.LewisN. E.PalssonB. O. (2011). Elimination of thermodynamically infeasible loops in steady-state metabolic models. Biophysical J.100 (3), 544–553. 10.1016/j.bpj.2010.12.3707

  • CrossRef
  • Google Scholar

• 12

  SchmidtA.KochanowskiK.VedelaarS.AhrneE.VolkmerB.CallipoL.et al (2016). The quantitative and condition-dependent Escherichia coli proteome. Nat. Biotechnol.34 (1), 104–110. 10.1038/nbt.3418

  • CrossRef
  • Google Scholar

• 13

  SchuetzR.KuepferL.SauerU. (2007). Systematic evaluation of objective functions for predicting intracellular fluxes in Escherichia coli. Mol. Syst. Biol.3, 119. 10.1038/msb4100162

  • CrossRef
  • Google Scholar

• 14

  WortelM. T.PetersH.HulshofJ.TeusinkB.BruggemanF. J. (2014). Metabolic states with maximal specific rate carry flux through an elementary flux mode. FEBS J.281 (6), 1547–1555. 10.1111/febs.12722

  • CrossRef
  • Google Scholar

• 15

  ZanghelliniJ.RuckerbauerD. E.HanschoM.JungreuthmayerC. (2013). Elementary flux modes in a nutshell: properties, calculation and applications. Biotechnol. J.8, 1009–1016. 10.1002/biot.201200269

  • CrossRef
  • Google Scholar