Law of Localization in Chemical Reaction Networks
In living cells, chemical reactions are connected by sharing their products and substrates, and form complex networks, e.g. metabolic pathways. Here we developed a theory to predict the sensitivity, i.e. the responses of concentrations and fluxes to perturbations of enzymes, from network structure alone. Responses turn out to exhibit two characteristic patterns, and . We present a general theorem connecting sensitivity with network topology that explains these characteristic patterns. Our results imply that network topology is an origin of biological robustness. Finally, we suggest a strategy to determine real networks from experimental measurements.
Cells have many chemical reactions, each of which is mediated by organic catalysts, enzymes. Reactions are not independent but connected and form complex networks. The dynamics of chemical concentrations are considered as the origin of physiological functions. However, dynamical behavior based on the network is not understood well.
One experimental approach to study such network systems is sensitivity analysis where the amount or activity of enzyme is perturbed and responses (concentrations of chemicals in the system) are measured Ishii () (see Fig. 1). However, the results of such experiments are very difficult to interpret, because theoretical criteria to evaluate the results of perturbations from network structures are not established.
There are other difficulties in understanding dynamical behaviors of reaction systems in biology. First, although a huge amount of information of reaction networks is available on databases KEGG (); Reactome (); BioCyc (), they provide no more than knowledge of identified reactions in biochemistry. It is possible that the information is incomplete, including many unidentified reactions or regulations. Second, in spite of the recent progress in biosciences, it is still difficult or almost impossible to determine quantitative details of the dynamics, such as functions for reaction rates, parameter values, or initial states.
In order to circumvent these difficulties, we introduce a mathematical method, named structural sensitivity analysis Mochizuki (); monomolecular (), to determine responses of chemical reaction systems to the perturbation of the enzyme amount or activity based only on network structure. From analyses we found that qualitative responses at a steady state are determined from information of network structure only. We also found that response patterns, e.g. distribution of nonzero responses of chemical concentrations in the network, exhibit two characteristic features, and depending on the structure of networks and position of perturbed reactions. Finally we found a general theorem connecting the network topology and the response patterns directly, and governing the characteristic patterns of responses. This theorem, which we call the law of localization, is not only theoretically important, but also practically useful for examining real biological systems. In the context of adaptation, there were some previous studies, which reported confined nonzero responses in specific systems Ni (); Steuer (); Drengstig (); He (). However, they did not find general laws of such response patterns, nor any topological conditions.
Ii Structural Sensitivity Analysis
We study concentration changes in a reaction system under perturbation of reaction rate parameters, assuming that the system is in a steady state Mochizuki (); monomolecular (). We label chemicals by and reactions by . A state of the system is specified by concentrations and obeys the following differential equations F1 (); F2 ()
Here, is called a stoichiometric matrix. is called a flux, which depends on metabolite concentrations and also on a reaction rate . We do not assume specific forms for , but assume that each is an increasing function of its substrate concentration;
Below, we abbreviate and emphasize nonzero as .
In this framework, enzyme knockdown of the -th reaction corresponds to changing the rate as (triangles in FIG. 2). By assuming steady state Palsson (); FBA1 (); FBA2 (); FBA3 (), the flux is expressed, in terms of a basis of , as , where is the dimension of the kernel and are coefficients depending on reaction rates. Under the -th knockdown, we have
The -th component of is also expanded as
From (3) (4), the response of steady state concentration (circles in FIG. 2) and flux (arrows in FIG. 2) to each perturbation is determined from network structure only Mochizuki (); monomolecular (). The result for each perturbation is simultaneously obtained through the following matrix computation:
where the matrix is given as
in a matrix notation. We call the inverse of as the sensitivity matrix . Note that , , , are column vectors with , and components respectively, and is an -by- matrix. We assume networks with throughout this paper, which guarantees the matrix is square, i.e. .
Comments are in order. First, our theory depends only on the structure of reaction networks. The network structure is reflected in the distribution of nonzero entries in the -matrix, which determines the qualitative responses. Second, as a generalization of our method, we account for regulations such as allosteric effects by relaxing (2) as
Then, regulations add additional in the -matrix, and the response is still determined through (5).
Iii Localization and Hierarchy
Let us see some results of structural sensitivity analysis.
: We consider a straight pathway, shown in FIG. 2 (Left). The -matrix and the sensitivity matrix are
The flux changes only when we perturb the top reaction 1 (the 1st column of ). The perturbation to reactions 2 or 3 changes only its substrate concentration (the 2nd, 3rd column of ).
: The second example shown in FIG. 2 (Right) consists of 6 reactions and 4 chemicals. The matrices A and are
Again, only the perturbation to the input rate, corresponding to the 1st column in (10), affect all chemicals and fluxes. Perturbations to reactions only decrease the concentrations of the substrates respectively. Knockdown of reaction decreases the concentrations along the cycle downward of the perturbation (see FIG. 2, and the 4th column of ). Knockdown of reaction 6 does not change the further downstream but does change in the cycle. Also, the signs of the responses are reversed (the 6th column of ).
: The third network in FIG. 3 (Left) includes 10 chemicals and 15 reactions. FIG. 3 (Right) shows nonzero response patterns of metabolites and inclusion relation between them. See Appendix for the -matrix and the sensitivity matrix. S
In general, response to perturbations in chemical reaction networks exhibits two characteristics, and . The means that the influence of the perturbations is confined in a finite region in a network. In other words, the naive intuition that a perturbation in an upper part of a reaction network influences all of the lower parts is incorrect. The implies that the nonzero response patterns under perturbations of different reaction rates exhibit inclusion relations among them.
Iv The Law of Localization
From the -matrix (6), we can generally prove a theorem, the law of localization, that determines the extent to which a perturbation influences in a network. For a given network, we consider a pair of a metabolite subset and a reaction subset satisfying the condition that includes all reactions influenced by metabolites in (see the condition (2’)). The choice of for a chosen is not unique in general. We call a subnetwork satisfying this condition “output-complete.” For such a subnetwork , we count the number of elements in , the number of elements in , and the number of the closed cycles that consist of the reaction subset . Then, we compute an index,
which is analogous to Euler characteristic and generally non-negative. The law of localization states that if for an output-complete subnetwork , then any perturbation of reactions in does not change the concentrations and the fluxes outside of , namely the perturbation effect is localized in itself. We call an output-complete subnetwork satisfying buffering structure.
Proof.– The theorem is proved from the distribution of nonzero entries of the -matrix. (i) Suppose a subnetwork is a buffering structure. Then by appropriately choosing a basis of the kernel of and the orderings of the indices of the -matrix, we can always rewrite the -matrix as
The lower left block vanishes because is output-complete. (ii) As explained already, the concentration change is proportional to , where is the minor matrix associated with the row of the -th reaction and the column of the -th metabolite. Then, for , follows because the upper left block in the minor , which was originally square in (12), is horizontally long.
: The network includes two buffering structures, and which are minimum buffering structures including only a single chemical and a single output reaction. The law of localization claims that the perturbation to reaction 2 in influence only the inside of , namely the concentration of A. Note the flux 2 in does not change in order to keep the outside of unchanged). We actually observed the predicted response in (8). Generally, a perturbation to a reaction which is a single output from a chemical influences the substrate concentration only.
: In addition to the 3 minimal buffering structures, , the network has two larger ones, (with ), (with ). is the minimum buffering structure including reaction 4. Then, the law of localization predicts that the nonzero response to perturbation of reaction 4 should be limited within , which is observed in the 4th column in (10). Similarly, the response to perturbation of reaction 6 is explained by .
: The network has 14 buffering structures, listed in Appendix. To examine the response hierarchy, we focus on the three buffering structures colored in FIG. 3; (with ), (with ), and (with ). Each of these three explains the response pattern under perturbation of reaction 5, 8, and 10 (or 13), respectively, and they satisfy an inclusion relation, . Accordingly, we can see from FIG. 3 (Right) that these response patterns satisfy an inclusion relation.
In this way, we understand all of the observed patterns from network topology by using the law of localization. In short, the first characteristic, , is explained from the existence of buffering structures. The second property, , is explained as the nest of the buffering structures.
Finally, as an application to real biological networks, we examine the carbon metabolism pathway of . The network is a major part of the energy acquisition process, and the basic structures are shared between bacteria and human beings. FIG. 4 shows the network Ishii (), including 28 metabolites and 46 reactions, and FIG . 5 shows the response hierarchy (see Appendix for the detail). Again, the response patterns exhibit the two characteristic features, and . The network has 17 buffering structures, and the existence and the nest of them explain the two characteristic features perfectly. We mention that some of the buffering structures, which are of course defined from network topology, surprisingly overlap biologically identified sub-circuits, the pentose phosphate pathway (yellow in FIG. 4, 5), the tricarboxylic acid cycle (blue) and the glycolysis (green). This correspondence may be understood from an evolutional point of view by considering the advantage of buffering structures.
V Discussions and Conclusions
Here we discuss the biological significances of buffering structures (and nest of them) in two different levels. The first discussion is on the physiological importance. A buffering structure prohibits influence of given perturbation from expanding to the outside, like a “firewall.” In other words, it is a substructure with robustness emerging from the network topology. The carbon metabolism network of E. coli possesses multiple nested firewalls (FIG. 5), and are expected to be robust to fluctuations of enzymes in it. We expect that such a topological characteristic of reaction networks could be the evolutionary origin of homeostasis in biological systems. A set of chemical reactions satisfying the condition of buffering structure by chance in evolutionarily early time would be positively selected as an advantageous circuit. We then expect that buffering structures in existing biological networks today might be generated and selected in such ways.
The second discussion is about practicality of the law of localization in experimental biology. Our knowledge of biochemical networks is considered incomplete: There might exist unidentified reactions or regulations. The condition for buffering structure depends on the local network structure only, which implies that we can study the sensitivity of the system only from local information on the network.
From this property, we can determine a “true” network by combining experiments as shown in FIG. 6. If a given perturbation (knock down or overexpression) to a predicted buffering structure, determined from network topology, influence outside of the buffering structure, then there must be inconsistency between the database information and the actual network. The mismatch must exist inside of the candidate structure, i.e. there must be unknown reactions or unknown regulations inside (or emanating from) the candidate subnetwork. By repeating theoretical predictions and experimental verifications, we can determine the “true” network from a partial network to the whole network in a step-by-step manner, i.e. from smaller to larger buffering structures. Our theory must promote the understanding of reaction networks in both the theoretical and experimental levels by directly connecting the network topology with behaviors of the systems.
Using a different method, Steuer et al. studied a mathematical criteria for “perfect adaptation,” where changing a rate constant in one part of the network does not affect steady-state concentrations or fluxes, which in fact, is a subpart of the phenomena we studied in this paper. There are at least three large differences: (i) We studied not only perfect adaptation, but also any qualitative responses (increase/decrease/invariant), (ii) While Steuer et al.’s method needs to examine a condition one by one for each pair of perturbation and chemicals, our method determines changes of all concentrations and fluxes by each perturbation of all reaction rates simultaneously via (5). (iii) We found and proved a general law which claims that the property of perfect adaptation emerges from local topology of network. Despite these differences, it would be interesting to explore relations between two mathematical theories.
This work was supported partly by the CREST, Japan Science and Technology Agency, and by iTHES research program RIKEN, by Grant-in-Aid for Scientific Research on Innovative Area, “Logics of Plant Development,” Grant number 25113005. We greatly appreciate Bernold Fiedler, Hiroshi Matano, and Hannes Stuke for their mathematical discussions. We also express our sincere thanks to Testuo Hatsuda, Michio Hiroshima, Yoh Iwasa, Sinya Kuroda, Masaki Matsumoto, Keiichi Nakayama, Madan Rao, and Yasushi Sako for their helpful discussions and comments.
- (1) Ishii, N., Nakahigashi, K., Baba, T., Robert, M., Soga, T., Kanai, A., Tomita, M. (2007). Multiple high-throughput analyses monitor the response of E. coli to perturbations. Science, 316(5824), 593-597.
- (2) Ogata, H., Goto, S., Sato, K., Fujibuchi, W., Bono, H., Kanehisa, M. (1999). KEGG: Kyoto encyclopedia of genes and genomes. Nucleic acids research, 27(1), 29-34.
- (3) Joshi-Tope, G., Gillespie, M., Vastrik, I., D’Eustachio, P., Schmidt, E., de Bono, B., Stein, L. (2005). Reactome: a knowledgebase of biological pathways. Nucleic acids research, 33(suppl 1), D428-D432.
- (4) Karp, P. D., Ouzounis, C. A., Moore-Kochlacs, C., Goldovsky, L., Kaipa, P., Ahrén, D., López-Bigas, N. (2005). Expansion of the BioCyc collection of pathway/genome databases to 160 genomes. Nucleic acids research, 33(19), 6083-6089.
- (5) Varma, A., Boesch, B. W., Palsson, B. O. (1993). Stoichiometric interpretation of Escherichia coli glucose catabolism under various oxygenation rates. Applied and environmental microbiology, 59(8), 2465-2473.
- (6) Edwards, J. S., Palsson, B. O. (2000). The Escherichia coli MG1655 in silico metabolic genotype: its definition, characteristics, and capabilities. Proceedings of the National Academy of Sciences, 97(10), 5528-5533.
- (7) Kauffman, K. J., Prakash, P., Edwards, J. S. (2003). Advances in flux balance analysis. Current opinion in biotechnology, 14(5), 491-496.
- (8) Orth, J. D., Thiele, I., Palsson, B. O. (2010). What is flux balance analysis?. Nature biotechnology, 28(3), 245-248.
- (9) Ni, Xiao Yu, Tormod Drengstig, and Peter Ruoff. The control of the controller: molecular mechanisms for robust perfect adaptation and temperature compensation. Biophysical journal 97.5 (2009): 1244-1253.
- (10) Steuer, Ralf, et al. Robust signal processing in living cells. PLoS Comput Biol 7.11 (2011): e1002218.
- (11) Drengstig, T., et al. A basic set of homeostatic controller motifs. Biophysical journal 103.9 (2012): 2000-2010.
- (12) He, Fei, Vincent Fromion, and Hans V. Westerhoff. Perfect robustness and adaptation of metabolic networks subject to metabolic and gene-expression regulation: marrying control engineering with metabolic control analysis. BMC systems biology 7.1 (2013): 131.
- (13) Mochizuki, A., Fiedler, B. (2015). Sensitivity of chemical reaction networks: A structural approach. 1. Examples and the carbon metabolic network. Journal of theoretical biology, 367, 189-202.
- (14) Fiedler B., Mochizuki A. (2015). Sensitivity of Chemical Reaction Networks: A Structural Approach. 2. Regular Monomolecular Systems. Math. Meth. Appl. Sci. 38, 3381-3600.
- (15) Craciun G., Feinberg M. (2006) Multiple equilibria in complex chemical reaction networks: The species-reactions graph. SIAM J. App. Math. 66:4, 1321-1338.
- (16) Feinberg M. (1995) The existence and uniqueness of steady states for a class of chemical reaction networks. Arch. Rational Mech. Anal. 132, 311-370.
Appendix A The analysis of Example 3
a.0.1 The matrix and sensitivity matrix for Example 3
The matrix for the network of Example 3 in the main text is given by
Here the row indices are the reactions , and the column indices are
where is a basis of the kernel space of the stoichiometric matrix . The vertical line separates the indices of the chemicals and those of the kernel vectors. By inverting , we obtain the sensitivity matrix . The result is , where are the diagonal matrices defined as
and is defined as
Here we have defined the factors () by
and the components in by
a.0.2 List of buffering structures
The network of Example 3 has the following fourteen buffering structures (and unions of them).
Appendix B E. coli central metabolism
b.0.1 List of reactions
1: Glucose + PEP G6P + PYR.
2: G6P F6P.
3: F6P G6P.
4: F6P F1,6P.
5: F1,6P G3P + DHAP.
6: DHAP G3P.
7: G3P 3PG.
8: 3PG PEP.
9: PEP 3PG.
10: PEP PYR.
11: PYR PEP.
12: PYR AcCoA + CO2.
13: G6P 6PG.
14: 6PG Ru5P + CO2.
15: Ru5P X5P.
16: Ru5P R5P.
17: X5P + R5P G3P + S7P.
18: G3P + S7P X5P + R5P.
19: G3P + S7P F6P + E4P.
20: F6P + E4P G3P + S7P.
21: X5P + E4P F6P + G3P.
22: F6P + G3P X5P + E4P.
23: AcCoA + CIT.
24: CIT ICT.
25: ICT 2KG + CO2.
26: 2-KG SUC + CO2.
27: SUC FUM.
28: FUM MAL.
29: MAL OAA.
30: OAA MAL.
31: PEP + CO2 OAA.
32: OAA PEP + CO2.
33: MAL PYR + CO2.
34: ICT SUC + Glyoxylate.
35: Glyoxylate + AcCoA MAL.
36: 6PG G3P + PYR.
37: AcCoA Acetate.
38: PYR Lactate.
39: AcCoA Ethanol.
40: R5P (output).
41: OAA (output).
42: CO2 (output).
43: (input) Glucose.
44: Acetate (output).
45: Lactate (output).
46: Ethanol (output).
b.0.2 List of buffering structures
The E. coli network exhibits the following 17 different buffering structures ().