Mathematical modeling of senescence in metabolic networks Ivanov, Oleksandr
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date: 2018
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Ivanov, O. (2018). Mathematical modeling of senescence in metabolic networks. University of Groningen.
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
senescence in metabolic networks
Printing: Ipskamp Printing
ISSN: 1570-1530
ISBN: 978-94-034-0519-3 printed version ISBN: 978-94-034-0518-6 electronic version
This Ph.D. study was carried out at the Groningen Institute for Evolu-tionary Life Sciences and at Johann Bernoulli Institute of Mathematics and Computer Science at the University of Groningen, the Nether-lands, and was financially supported by the Netherlands Organisation for Scientific Research (NWO).
Copyright c by Oleksandr Ivanov. All rights reserved. No part of this thesis may be reproduced, stored in a retrieval system or transmitted in any form or by any means without the prior written permission of the author.
senescence in metabolic networks
PhD thesis
to obtain the degree of PhD at the University of Groningen
on the authority of the Rector Magnificus, Prof. E. Sterken
and in accordance with
the decision by the College of Deans. This thesis will be defended in public on
Monday 5 March 2018 at 14.30 hours
by
Oleksandr Ivanov
born on 17 September 1987 in Berezne, Ukraine
Assessment commitee Prof. dr. B. M. Bakker Prof. dr. B. Jayawardhana Prof. dr. J. Molenaar
1 Introduction and thesis overview 1
1.1 General models of metabolic networks . . . 3
1.2 Network topology . . . 4
1.3 Network function (optimal flux) . . . 6
1.4 Metabolic control analysis . . . 8
1.5 Properties of steady states and dynamics of metabolite concentrations . . . 9
1.6 Thesis overview . . . 11
2 A guide to chemical reaction network theory 15 2.1 Introduction . . . 16
2.2 Chemical reaction networks . . . 17
2.2.1 Structure . . . 19
2.2.2 Kinetics . . . 20
2.3 Dynamical properties of CRNs . . . 22
2.4 General kinetics . . . 23
2.4.1 Stoichiometric and Jacobian matrices . . . 23
2.4.2 Species-reaction graph . . . 24
2.5 Mass action kinetics . . . 26
2.5.1 Stoichiometric and Jacobian matrices . . . 29
2.5.2 Balancing . . . 29
2.5.3 Deficiency . . . 30
2.6 Concluding remarks . . . 33
3 Steady states and their stability in metabolic networks without regulation 35 3.1 Introduction . . . 37
3.2 Metabolic networks . . . 39
3.2.1 Topology and stoichiometry . . . 40
3.2.2 Reaction kinetics . . . 42 3.2.3 Linear approximation . . . 44 3.3 Steady states . . . 46 3.4 Single-substrate-single-product networks . . . 49 3.5 Multiple-substrates-multiple-product networks . . . 61 3.5.1 Simple stoichiometry . . . 61 3.5.2 General stoichiometry . . . 66 3.6 Discussion . . . 75 3.7 Appendix A . . . 78 3.8 Appendix B . . . 81 3.9 Appendix C . . . 83 3.10 Appendix D . . . 93 3.11 Appendix E . . . 96
4 Potential effects of senescence on metabolic networks with regulation 99 4.1 Introduction . . . 100
4.2 Metabolic systems without regulation . . . 101
4.3 Alternative steady states . . . 104
4.4 Oscillations . . . 106
4.5 Oscillations and chaos . . . 108
5 Evolution of metabolic control in the absence and pres-ence of senescpres-ence 113 5.1 Introduction . . . 114
5.2 The model . . . 116
5.2.1 Overview of the model structure . . . 116
5.2.2 Flux regulation . . . 117
5.2.3 Evolution . . . 119
5.3 Results . . . 120
5.3.1 Fitness landscapes . . . 120
5.3.2 Evolutionary simulations . . . 127
5.3.3 The influence of senescence . . . 134
5.3.4 Degradation of baseline enzyme activity . . . . 134 6
5.3.5 Degradation of regulatory coefficients . . . 137
5.4 Appendix A . . . 141
5.5 Appendix B . . . 144
5.6 Appendix C . . . 146
6 Summarizing discussion 149 6.1 Metabolism and senescence . . . 150
6.2 My contribution to the discussion . . . 153
6.3 Robustness of metabolic networks . . . 155
6.4 Criticism . . . 156 6.5 Future perspectives . . . 157 Bibliography 160 Samenvatting 183 Acknowledgements 187 7
Introduction and thesis overview
T
his thesis began as a project out of the Systems Biology Centre for Metabolism and Ageing at the University of Groningen. The center aims to understand the connections between metabolism and senescence, particularly, how metabolic networks deteriorate with age and the effects of that deterioration on the function of the organism as a whole. The aim of my project was to provide a theoretical background for understanding senescence in metabolic networks.Senescence can be defined as a decline in the functioning of living systems. Organismal senescence is characterized by a declining ability to respond to stress, increased homeostatic imbalance, and increased risk of disease. A wealth of information has been collected about senes-cence during decades of investigation. The Encyclopedia of senessenes-cence, which contains nearly 2000 pages, only partially covers the information known about senescence [1]. By the year 1990, there were over 300 the-ories of senescence [2], and this number continues to grow. Perhaps, there is no other area in the life sciences to which such a large number of theories is dedicated.
Many factors have been reported to influence senescence. Among them, the most prominent are free radicals [3, 4, 5, 6, 7, 8], inflam-mation [9], replication stress [10], various toxins [11] and viruses [12], temperature [13, 14], telomere shortening [15, 16], and regulation of gene expression [11]. These factors lead to the degradation of DNA, proteins, lipids, and membranes, resulting in senescence.
During the last century, senescence research has been dominated by a reductionist approach [17]. Most effort was spent on identification
and characterization of individual genes and proteins that modulate senescence [18]. However, cell components carry out their functions within complex networks, where a single component can affect a wide range of other components [19]. It has become apparent that senes-cence is a systemic property and involves various processes that often operate in parallel. As a result, systems biology approaches (i.e., ge-nomics, transcriptomics, proteomics, metabolomics) have begun to be used in senescence research [20, 21, 22, 23, 24]. Therefore, to inves-tigate senescence in a systematic manner, research should focus on understanding senescence processes at the network level. An effort to approach senescence research in this manner is already underway. Specifically, genetic determinants of longevity in Saccharomyces cere-visiae have been identified using network analyses [25]. Intracellular [26] and matrix [19] protein-protein interaction networks related to senescence have also been investigated. Any change in gene regulation, protein function, and gene regulatory or protein-protein interactions ultimately reflects on metabolism and its regulation. Therefore, the metabolic network is the ultimate biological network that experiences senescence [24].
Generally, there are two ways to investigate senescence in metabolic networks theoretically. The first approach is to build detailed models of specific pathways and investigate how deterioration affects the func-tions performed by these pathways. Benefits of this approach are that the structures of many metabolic pathways are well known, and it re-quires relatively small-scale modeling. The downside to this approach is that even with well-known pathways, the specific kinetics of each reaction are not always known and in most cases, few parameters are known. Therefore, the accuracy of detailed modeling critically depends on assumptions made about reaction kinetics and parameters.
The second approach is to build general models of metabolic net-works in an attempt to understand their general properties and how they are affected by deterioration. In this case, the aim is to derive properties of metabolic pathways or networks without detailed knowl-edge of reaction kinetics and parameters.
For my thesis, I chose the second approach, building general mod-els, to investigate senescence in metabolic networks. The fundamental
properties I considered were the number and stability of steady states, the dynamics of metabolite concentrations, and the roles of metabolic regulation and control. The choice to build general models rather than models of specific pathways was motivated by the idea that senescence in metabolic networks is related to their properties. This idea does not exclude the possibility that senescence is related to individual prop-erties of metabolic pathways. In other words, I did not reject the idea that detailed modeling of specific pathways is useful for modeling senescence. The advantage to using general models is that discoveries may provide a bird’s-eye perspective on metabolism and outline new and interesting directions in the research of metabolic networks and senescence.
General properties of metabolic networks are closely connected to network robustness. Robustness is a property that allows a system to maintain its functions despite external and internal perturbations [27] and is considered to be one of the fundamental features of complex evolved systems [27, 28, 29]. A counterpart to robustness, fragility, is related to senescence. In my thesis I will show that robustness of some metabolic network properties can lead to fragility and, as a result, to potential senescence in other network properties.
1.1
General models of metabolic networks
Different modeling approaches can be used to overcome a lack of de-tailed knowledge about the metabolic reactions, while still drawing general conclusions about the metabolic network properties and ro-bustness.
The first approach considers only the topology of metabolic net-works. There are various ways to construct a graph from the structure of a metabolic network, which will be discussed in the section Network topology. Network topology alone contains valuable information about the organization of networks and their robustness.
A second approach is based on the assumption that metabolism is at a steady state and metabolic networks optimize a combination of specific fluxes to achieve optimal functionality, such as optimal growth
(in bacteria and yeast) or ATP production. This approach is the foun-dation of the flux-balance analysis [30].
Another approach considers the sensitivity of certain network prop-erties (e.g., fluxes or metabolite concentrations) to changes in cer-tain parameters (e.g., enzyme activity). This approach is known as metabolic control analysis [31].
The approaches discussed above exclude metabolite concentration dynamics, consideration of which requires involvement of reaction ki-netics. As I will demonstrate in Chapters 2-4 of my thesis, it is possible to draw general conclusions about the dynamics of metabolic networks using general assumptions about reaction kinetics and without knowl-edge of specific parameter values.
I discuss these approaches below in detail.
1.2
Network topology
Topological approaches use a graph perspective to model metabolic networks.
One topological approach considers metabolites as nodes. Reac-tions represent connecReac-tions (links) between those nodes. Two nodes are connected if they participate in the same reaction. The graph of a metabolic network is built upon such nodes and links. Jeong et al. considered the topology of metabolic networks using this approach; they used genome data from 43 sequenced organisms to reconstruct the metabolic networks of those organisms [32]. Jeong et al. proposed that all 43 metabolic networks were scale-free [32], meaning that the network degree distribution (i.e., the number of links per node) fol-lowed a power-law distribution, Ppkq ky. Here, Ppkq is the fraction of nodes in the network having k links and y is a parameter. Later scale-free structure of reconstructed metabolic networks of 80 organ-isms was demonstrated by Ma and Zeng [33]. This result can be ap-plied to assess the robustness of metabolic networks. It is known that in a scale-free network, the network diameter, defined as the average length of the shortest path between any two nodes in the network, does not increase substantially after random removal of nodes [34].
How-ever, the targeted removal of highly-connected nodes (hubs) leads to a rapid increase in network diameter [34]. That is, the network diameter of scale-free networks is robust against random attacks but vulnera-ble to targeted attacks. This result was also shown for reconstructed metabolic networks [32]. In addition, random removal of nodes from scale-free networks allows for the preservation of the giant connected component (the largest set of connected nodes), while the number of smaller connected components increases [34]. In contrast, targeted re-moval of highly-connected nodes leads to disintegration of the network to its components [34]. Complementary to this, Ferrarini et al. [35] showed that restoration of hubs in a scale-free network restores 35% of the lost function, while restoration of weakly-connected nodes did not have a significant effect. The idea that metabolic networks have a scale-free structure was disputed by several authors [36, 37].
In addition to the scale-free, a bow-tie structure of metabolic net-works was proposed [38]. Bow-tie structure is represented by a con-servative core (the knot of the bow-tie) that consists of the essential and highly-connected metabolites of all species (e.g., ATP, NADP, glu-cose, nucleotides, amino acids and others). Nutrients fan into the core and building blocks, such as proteins, DNA, fatty acids, and sugars fan out [38]. Interestingly, the bow-tie itself consists of other bow-ties. For example, the core metabolism consists of precursors from which other precursors are made (nucleotides, amino acids, sugars and oth-ers) through anabolism. A bow-tie structure is proposed to promote metabolic robustness because the high connectivity of the core allows for alternative routes if nodes are damaged and because of the ability to change building blocks within the core [38]. In an analysis of 65 re-constructed metabolic networks, Ma and Zeng supported the bow-tie structure of metabolic networks [39]. Bow-tie and scale-free structure of metabolic networks are not mutually exclusive [39].
Another network topology approach considers metabolic networks as bipartite graphs. In a bipartite graph, there are two types of nodes: metabolites and reactions. A metabolite node is connected to a reac-tion node if the metabolite participates in the corresponding reacreac-tion. Lemke et al. constructed a bipartite graph from the annotated genome sequence of Escherichia coli [40] and used it to assess network
robust-ness. They defined network damage as the number of metabolites that cannot be produced as a result of removing the corresponding reaction. Lemke et al. showed that 91% of enzymes cause little or no damage when removed, whereas 9% can cause substantial damage [40]. Additionally, a study of young (2 years old) and old (13 years old) marmosets showed that the connectivity of a metabolic network decreases with age [24], which supports the modeling of the mechanism of damage in the study by Lemke et al.
The studies above indicate the crucial role of hubs in network ro-bustness. Some authors also argue that weakly-connected nodes may be just as important to robustness and that failure of these nodes may lead to the loss of network integrity, thus playing a role in senescence [41].
Many other topological features of metabolic networks have been in-vestigated. In particular, it was proposed (and disputed) that metabolic networks are small-world [32, 42], highly clustered (modular) [43], dis-assortative (i.e., nodes with high connectivity tend to attach to nodes with low connectivity) [43], hierarchical [44], and self-similar [45] with conserved hubs (i.e., most connected nodes) [32] and highly conserved network diameters between species [32].
According to chemical reaction network theory (CRNT), reaction complexes (the set of substrates or products of a given reaction) are considered as nodes of metabolic networks and reactions linking those nodes. Chemical reaction network theory will be discussed in Chapters 2 and 3 in detail.
1.3
Network function (optimal flux)
Although evaluating the robustness of a metabolic network using purely topological considerations can provide insight into the potential mech-anisms of function and degradation, it does not address metabolism dynamics. To address this topic, flux-balance analysis investigates the rate of change of metabolite concentrations (fluxes) at steady state, during which concentrations of metabolites and fluxes are constant [30]. That is, if x is a vector of metabolite concentrations, then the
rate of change of x at steady state is zero: dx
dt 0. (1.1)
In flux-balance analysis, it is assumed that the goal of a metabolic network is to optimize (in most cases to maximize) a specific function z (called the objective function), subject to flux constraints. The ob-jective function is defined as a weighted sum of specific fluxes (reaction rates):
z c1v1 c2v2 . . . cnvn. (1.2)
Here, c1, c2, . . . , cn are constants and v1, v2, . . . , vn are fluxes that
are used to assess the objective function.
Examples of common objective functions are the rate of ATP pro-duction and the sum of fluxes of reactions that produce precursor metabolites for cell growth [30].
The objective function is optimized given the following constraints on fluxes:
vi,min ¤ vi ¤ vi,max, i 1, ..., n. (1.3)
A flux balance analysis was used to assess the robustness of metabolic networks in various organisms. In this case, robustness was defined as the ability of a metabolic network to achieve the objective function un-der knockout of metabolic reactions or change of flux constraints. One study, which analyzed the metabolic capabilities of E. coli, determined that inhibition of the majority of reactions would not substantially affect growth. By contrast, there are a small number of essential re-actions (14% of all rere-actions) for which inhibition by gene knockout would lead to a substantial decrease in the growth objective function [46]. Significantly reducing flux through many of these essential re-actions also showed little effect on growth as the objective function [47]. Similar results were obtained for Haemophilus influenzae [48]. It was shown for this organism that metabolic networks achieve ro-bustness through redistribution of fluxes through undamaged parts of the network [46, 47]. Another study of Saccharomyces cerevisiae using
flux-balance analysis showed that robustness of cell growth is related both to hub nodes (i.e., those with a large number of links to other metabolites) and metabolites with a low number of links that represent essential steps in certain pathways [49]. These findings are in qualita-tive agreement with the results obtained using topological approaches to identify the robustness of metabolic networks.
1.4
Metabolic control analysis
Metabolic control analysis investigates how network properties (e.g., flux or metabolite concentrations) depend on the properties of network components (e.g., the activities of individual enzymes). In other words, metabolic control analysis studies the sensitivity of network properties to the properties of network components. The sensitivity is described formally by the concept of control coefficients. For example, sensitivity of flux in a pathway to activity of the enzymes that constitute the pathway is described by the flux control coefficient:
CiF BF {F BEi{Ei
. (1.4)
Here, F is the total flux in the pathway of interest, and Ei is the
activity of an enzyme i [31]. If the control coefficient of a certain enzyme is equal to zero, then the enzyme has no control over the pathway; if the control coefficient is close to one, then the enzyme has a strong control over the pathway. Several general properties of metabolic control have been discovered. One, the summation theorem, states that in any metabolic pathway the sum of all control coefficient is always equal to one [31]. This theorem has significant consequences for understanding how systemic properties in metabolic pathways (such as flux) are controlled. For example, in linear metabolic pathways if one control coefficient increases, the others must decrease. Therefore, each control coefficient has a system-wide influence over the activity of individual enzymes.
Historically, it was thought that there is one enzyme (with activity much lower than that of other enzymes) controlling the flux through a pathway. This enzyme was called the rate-limiting enzyme [50, 51].
According to the summation theorem, the rate-limiting enzyme has a flux control coefficient close to one, whereas the other enzymes in the system have control coefficients close to zero. However, experi-mental evidence shows that true rate-limiting enzymes are rare [31]. The modern view is that control over a pathway is distributed among all enzymes [52, 53]. Thus, the question about pattern of control over the pathway (i.e., the distribution of control coefficient values) remains open. Chapter 5 will address this question and discuss how the distri-bution of control coefficients contributes to the robustness of metabolic pathways.
1.5
Properties of steady states and dynamics of
metabolite concentrations
Steady-states and dynamics of metabolite concentrations are another important aspects of metabolic network function. Important prop-erties of steady states are the number of steady states (i.e., unique steady states, multiple separated steady states, manifolds of steady states) and their stability (locally stable, globally stable, unstable). The location of a steady state in the space of metabolite concentrations determines the state of metabolism. For example, proper concentra-tions and production rates of metabolites, such as ATP, AMP, NAD, fructose-6-phosphate, are crucial for proper cell functioning. A change in the concentrations of certain metabolites (e.g., neurotransmitters, fatty acids, glucose) may lead to various diseases and senescence. The way a metabolic system copes with deviations in metabolite concen-trations after a perturbation of the system is of crucial importance and is described by the system’s dynamics and stability of steady states.
Modeling steady state properties and metabolite concentration dy-namics requires consideration of the full complexity of metabolic net-works. That is, a large number of metabolites (represented in models by dynamic variables) translates to a large number of differential equa-tions. Moreover, complex network connectivity and the presence of regulation contribute to complexity of differential equations. Adding to this complexity, the form of the reaction rates and the values of
One way of managing metabolism complexity is to investigate small subnetworks for which data is available regarding reaction rates and parameters. To this end, much information is known about the dy-namical behavior of specific metabolic pathways (e.g., glycolysis, TCA cycle, fatty acids metabolism) [54, 55, 56, 57]. Models of these specific pathways are able to explain various phenomena, such as glycolytic oscillations [58, 59].
Moreover, full-scale dynamic models of entire metabolic networks have been developed (e.g. for the human erythrocyte) [60, 61, 62]. However, these models shed little light on general properties of steady-states and dynamics of metabolite concentrations.
Generally, dynamics of metabolite concentrations can be described by the following equation:
dx
dt Svpxq. (1.5)
Here, x is the vector of metabolite concentrations, S is the stoichio-metric matrix, and vpxq is the vector of reaction rates. In the stoichio-metric matrix, rows correspond to different metabolites and columns correspond to different reactions. Each entry in the stoichiometric ma-trix represents the number of metabolite molecules that participate in a certain reaction. Equation (1.5) is fundamental for understanding metabolic network dynamics. In particular, equation (1.5) splits the representation of metabolic network dynamics into two parts: a struc-tural part, represented by the stoichiometric matrix, S and a kinetics part, represented by the vector of reaction rates, vpxq.
An approach that relies on general assumptions of network struc-ture and reaction kinetics is a productive way to investigate steady-state properties and metabolic dynamics in the absence of detailed knowledge about reactions. This direction of research is known as the chemical reaction network theory [63]. Important results of the chem-ical reaction network theory are reviewed in Chapter 2. A current criticism of this theory is that many results are based on restrictive assumptions in relation to biological realism, e.g. mass action kinetics or the necessity of inflows and outflows of all metabolites. In Chapter
3, I take this concept one step further by deriving various steady-state properties based on very general and realistic assumptions regarding reaction kinetics and network structure.
1.6
Thesis overview
Despite the challenges of modeling metabolism, a number of solid the-oretical approaches have been developed that allow the investigation of general metabolic network properties. Using the approaches outlined above, progress have been achieved in understanding the topological structure of metabolic networks, network function, control of network fluxes, steady-state stability, metabolite concentration dynamics, and others. However, the results generated by each of these approaches and even the theories themselves are currently topics of active debate. In this thesis, I took a step toward understanding the general prop-erties of metabolic networks with an emphasis on robustness and senes-cence. I started by using the chemical reaction network theory to in-vestigate steady state number, stability, and dynamics. I used this approach to identify ’spots of vulnerability’ in metabolic networks. A spot of vulnerability is the part of a metabolic network in which steady-state properties or dynamical behavior change as a result of parameter changes. For simplicity, I started by considering steady-state proper-ties and dynamics in metabolic networks without regulation (Chapters 2 and 3).
I then considered senescence in metabolic networks with regulation (Chapter 4). In Chapter 5, I used metabolic control analysis to in-vestigate the general principles of regulation and control in metabolic networks in the context of senescence.
The most important insight I gained through this work was that properties of metabolic networks (such as regulation) that provide ro-bustness to one aspect of network functioning often lead to ity in other aspects of functioning. As a consequence, such vulnerabil-ities may lead to senescence. I discuss this in Chapter 6 in detail.
In Chapter 2, I review important developments and conclusions of chemical reaction network theory. Using very general assumptions
theory is able to draw general conclusions about steady-state prop-erties and metabolite concentration dynamics. A main result of this theory is that in most cases, metabolic networks without regulation are expected to have a unique global, asymptomatic steady-state stability or a globally stable steady-state manifold. As a result, the dynamical behavior of such systems is very simple, namely, a monotonic conver-gence to a steady state.
Despite its generality, a problem with the chemical reaction network theory is that some of its assumptions are not biologically realistic. For example, most of the results of this theory are based on the assumption of mass action, which rarely occurs in biological systems. Another unrealistic assumption is that the metabolic system is considered as a continuous flow stirred-tank reactor, meaning that all metabolites have inflow and outflow reactions, while the total mass is conserved. A more realistic assumption would be that certain metabolites (i.e., precursors) flow into the network and certain metabolites (products and waste) flow out.
In view of this, in Chapter 3 I contributed to the chemical reac-tion network theory by considering more realistic assumpreac-tions regard-ing network functionregard-ing. First I investigated steady-state properties without making any assumptions about reaction kinetics. Then, to investigate steady-state properties and metabolic dynamics in more detail, I considered two-way weakly monotonic reaction kinetics. That is, kinetics wherein the reaction rate increases with an increase in sub-strate concentration and does not increase with an increase of product concentration. Typical examples of reaction kinetics (e.g., mass action, Michaelis-Menten, and Hill kinetics) are two-way weakly monotonic. Importantly, monotonicity of kinetics excludes metabolic regulation. I also assumed that only some metabolites flow in and out of the system. Despite the generality of these assumptions, I drew conclusions about steady-state properties and metabolic dynamics.
Even when using realistic assumptions, models of metabolic net-works without regulation still demonstrate very simple dynamical be-havior, namely monotonic convergence to the set of steady states. In such metabolic networks, senescence that acts on enzymes and results
in the change of model parameters can only lead to changes in the location of steady states, not the dynamical behavior.
As a logical consequence, in Chapter 4 I modeled how senescence can act on metabolic systems with regulation. I considered simple examples of metabolic systems with regulation and introduced senes-cence. It is known that enzyme activities change over a lifetime, leading to changes in their reaction rates. I modeled senescence as a gradual change in parameters of reaction rates with time. The first general conclusion is that with regulation, even very simple metabolic systems are capable of rich dynamical behavior (e.g., oscillations, bistability, chaos) compared to unregulated metabolic systems. The second gen-eral conclusion is that regulated metabolic systems can dramatically change their dynamical behavior (undergoing bifurcations), as a result of gradual changes in parameters. In Chapter 4 I also discuss the im-plications of such dynamical behavior to the robustness of metabolic networks and senescence.
In Chapter 5 I considered the interplay of senescence with an-other aspect of metabolic network functioning, metabolic control. I developed an individual-based evolutionary model to study the opti-mal distribution of enzyme activities (represented by control coeffi-cients) within linear pathways in the absence and presence of senes-cence. Specifically, I considered a population of individuals with linear metabolic pathways. I assumed steady state metabolism and calcu-lated the overall flux in the pathway. I also assumed that there is an optimal flux that the pathway needs to achieve in order to attain maximal fitness. At certain time points, the flux changes randomly between high and low values. Metabolic regulation changes enzyme activity to optimize the flux. The time that it takes for regulation to achieve optimal flux depends on baseline enzyme activity and the strength of the regulatory processes. I evolved baseline enzyme activi-ties and regulation strength in order to study the optimal distribution of enzyme activities and the interplay between metabolic regulation and control. I modeled senescence as the degradation of enzyme activ-ities and regulation strength over an individuals lifetime to investigate how it affects the evolution of regulation and control.
pre-insight from my research is that the presence of metabolic regulation may be an important reason for metabolic network senescence. A sec-ond insight is that senescence in metabolic networks, and perhaps in whole organisms, can be the result of a trade-off between different as-pects of metabolic robustness. In brief, metabolic regulation provides robustness against deviations from optimal metabolite concentrations and fluxes on daily basis. However, it provides the background for vulnerability of metabolic networks to changes in reaction rate param-eters, which may lead to complex dynamical behavior and undesired metabolic states. These insights are presented in the light of the main theories that connect senescence and metabolism and modern views on robustness in the life sciences, which are discussed in Chapter 6 in detail.
A guide to chemical reaction network
theory
Abstract
Biochemical networks form the basis of life. Research ef-forts to model biochemical network function are hampered by the complexity of such networks. Complexity arises from the presence of a large number of interactive com-ponents, complex network topology, and unknown reac-tion rates form and parameters. Components (chemical reactants) and topology of biochemical networks is often known, whereas the specific forms of reaction rates and pa-rameters are often unknown. Although various theoretical approaches explore the available knowledge about network components and topology, they do not model the dynamics of those network components. Chemical reaction network theory (CRNT) is a theoretical approach that combines knowledge of the biochemical network structure with mild assumptions about reaction kinetics to investigate chemi-cal reaction dynamics in the absence of known parameters. Several technical reviews on CRNT are available, primarily for mathematically-oriented researchers. In this chapter, I provide a less formal guide to CRNT. Specifically, I present results on a number of steady states and their stability in chemical reaction networks.
2.1
Introduction
B
iochemical networks, which comprise metabolic, cell signaling, and regulatory networks, are crucial for understanding the func-tion of cells and often whole organisms. In the last two decades, there has been an intense research effort aimed towards understanding bio-chemical networks. This research must consider the complexity of such networks [33, 64, 65], which arises due to a large number of interact-ing components and intricate connections within and between different parts of the network. Many detailed models have been introduced to explain the function of specific biochemical networks and to assist ex-perimental research [64, 65, 66, 67]. However, it is often not clear whether a network model adequately represents the network because the kinetics of biochemical reactions is often unknown. In cases where kinetic laws are well-established, parameters are typically not known. In contrast, due to the rapid development of molecular and systems biology methods and genome-scale reconstructions [68], the topologi-cal structure of many biochemitopologi-cal networks is known. The question therefore arises, whether and to what extent such structural informa-tion is sufficient to provide insights into the dynamical properties of biochemical networks.Various theoretical approaches that rely on knowledge about the structure of biochemical networks have been developed to draw con-clusions about network properties without detailed knowledge about reaction kinetics and parameters. These ’structural’ approaches in-clude constraint-based modeling [69], structural kinetic modeling [70], applications of graph theory [32, 71, 72], and Boolean and rule-based models [66]. The disadvantage of most of the structural methods is that they do not explicitly include network dynamics. CRNT is an approach that uses the structure of chemical (and biochemical) net-works to investigate the dynamical behavior of the concentrations of their components.
Several technical reviews on developments in CRNT are available [63, 73, 74]. However, they primarily concentrate on mass action ki-netics and are oriented to a mathematical audience. The aim of this review is to provide a less technical introduction to CRNT and, in
particular, to highlight recent developments in this field. Therefore, in this review I address only some features of the dynamics of chemical reaction networks (CRNs), namely, the uniqueness of steady states, multistationarity (the presence of multiple steady states), and the sta-bility of steady states.
This chapter is structured as follows: in sections 1 and 2, I provide the most essential definitions related to CRNT; in section 3, I present recent results of CRNs on a number of steady states and their stability using general assumptions about kinetics; in section 4, I present re-sults assuming mass action kinetics; finally, in the concluding remarks section, I summarize and discuss these results in detail.
2.2
Chemical reaction networks
A chemical reaction network (CRN ) is a set of chemical species and the chemical reactions in which they participate. I begin this guide to CRNT with an illustrative example of a CRN, which will be used as a running example. Consider the following set of chemical reac-tions, which represent the phosphorylation-dephosphorylation process (an example of a futile cycle) [75] and often participate in signaling cascades, bacterial two-component systems, and other biological net-works: E S v1 ÐÑ ES v2 ÝÑ E P F P ÐÑ F Pv3 ÝÑ Fv4 S. (2.1)
Here, E, S, P, F, ES, and F P are chemical species. S is a sub-strate, P is a product, E and F are enzymes, and ES and F P are enzyme-substrate complexes. The concentrations of these species are denoted as xE, xS, xP, xF, xES, and xF P, and the reaction rates as
v1, v2, v3, v4. Reaction rates typically depend on the concentrations of
metabolites.
The concentration dynamics of these chemical species can be de-scribed by the following set of differential equations:
9 xE v2 v1 9 xS v4 v1 9 xF v4 v3 9 xP v2 v3 9 xES v1 v2 9 xF P v3 v4. (2.2)
In matrix form, this equation is written in the following way: 9 xE 9 xS 9 xF 9 xP 9 xES 9 xF P 1 1 0 0 1 0 0 1 0 0 1 1 0 1 1 0 1 1 0 0 0 0 1 1 v1 v2 v3 v4 . (2.3)
Here, 9xi dxi{dt. The stoichiometric matrix encapsulates the
structure of the CRN. The vector of chemical species is denoted as x pxE, xS, xF, xP, xES, xEPq and the vector of reaction rates is
de-noted as v pv1, v2, v3, v4q. The matrix that premultiplies the vector
v is denoted as S. A compact version of the dynamics of reactant con-centrations within the network (2.1), and in fact within any network [76], can then be represented as:
9x Svpxq. (2.4)
For a CRN with m chemical species and n reactions in the equation (2.4), x represents the m-dimensional vector of metabolite concentra-tions, S represents the mn stoichiometric matrix, and vpxq represents an n-dimensional vector of reaction rates.
Importantly, equation (2.4) allows chemical reaction dynamics to be divided into a structural part, represented by the stoichiometric matrix, and a kinetic part, represented by the vector of reaction rates vpxq. As will be demonstrated, the structure of a network (represented by the stoichiometric matrix S ) determines its dynamical properties using general assumptions regarding reaction kinetics.
2.2.1
Structure
In the stoichiometric matrix, rows correspond to different species and columns correspond to different reactions. The entry Sij represents the
number of molecules of chemical species i that participate in reaction j. If Sij 0, then i is a substrate of the reaction j ; if Sij ¡ 0, then
i is a product of the reaction. If Sij 0, then i does not participate
in the reaction j. If, in every reaction, there is only one molecule that participates, then the entries of the stoichiometric matrix are either 1,1, or 0. In the literature such stoichiometry is called binary [77]. The network (2.1) has binary stoichiometry.
The set of substrates or the set of products for a given reaction is called the reaction complex. For example, in the CRN (2.1), the reaction complexes are E S, ES, E P, F P, F P, and F S. These reaction complexes will be denoted as C1, . . . , C6. In terms of reaction
complexes, network (2.1) is of the form: C1 v1 ÐÑ C2 v2 ÝÑ C3 C4 v3 ÐÑ C5 v4 ÝÑ C6. (2.5) If every complex in a CRN contains precisely one chemical species, then such a network is called a single-substrate-single-product network. Otherwise, the network is called a multiple-substrate-multiple-product network.
CRNs are often represented as graphs, for which chemical species or reaction complexes are considered nodes and chemical reactions are considered links that connect those nodes.
A path between nodes i and j in a network is a sequence of links that leads from node i to node j.
A complex Ci is said to be connected to complex Cj if there exists
a path between them. Note that path and connectedness do not imply directionality.
A linkage class is a maximal set of connected reaction complexes. In network (2.5), there are two linkage classes: L1 tC1, C2, C3u and L2
tC4, C5, C6u.
Another aspect of a CRN is reversibility of reactions. A reversible chemical reaction can proceed in forward and backward directions,
The reversible network is one that consists of reversible reactions only. A CRN is weakly reversible if the existence of a path from complex Ci to complex Cj implies the existence of a path from complex Cj to
complex Ci. For example, consider the following CRN:
C1 Ø C2
Ò Ö Ò C5 Ø C6.
C3 Ñ C4
(2.6)
Although network (2.6) is not reversible and every complex is not connected to every other complex, it is considered weakly reversible.
Another aspect of network structure is inflow and outflow rules of the reactants, which play an important role in determining network properties. The simplest case is when the network is closed, meaning there is no inflow or outflow. In contrast, an open network allows inflow and outflow of chemical species. The most natural models allow the inflow and outflow of some reactants [78, 79, 80, 81]. Inflow rates are often assumed to be constant. A Continuous flow stirred-tank reactor (CFSTR) is a type of open network in which there is a constant inflow rate and a linear degradation (outflow) rate for all chemical species. Moreover, the degradation rate is the same for all chemical species and the sum of inflows and outflows is equal to 0 (i.e., total mass is conserved) [73]. A generalization of a CFSTR is a model in which some species are free to enter and leave the system, whereas others are entrapped [82, 83].
2.2.2
Kinetics
The kinetics of a CRN corresponds to the relationship between the rate of the reaction and the concentrations of the chemical species, which are described by the vector of reaction rates vpxq in equation (2.4). The simplest form of kinetics is mass action kinetics. For an irreversible reaction j with mass action kinetics, the reaction rate is:
vjpxq k
¹
s
xSsj
Here, xs are concentrations of reaction substrates and Ssj are the
entries of the stoichiometric matrix that define the number of molecules of substrate s that participate in reaction j.
For a reversible reaction j with mass action kinetics, the reaction rate is: vjpxq kf ¹ s xSsj s kr ¹ p xSpj p . (2.8)
Here, kf and kr are constants, xp are concentrations of reaction
products, and Spj are the entries of the stoichiometric matrix that
define the number of molecules of product p that participate in reaction j.
In biochemistry, Michaelis-Menten and Hill kinetics are often as-sumed.
The reaction rate with Michaelis-Menten kinetics is: vpxq vmaxxs
xs Km
. (2.9)
.
Here, vmax is the maximal reaction rate and Km is the Michaelis
constant. The reaction rate with Hill kinetics is: vpxq x n s xn s Kn . (2.10) .
Here, n is the Hill coefficient and K is a constant.
These types of kinetics are similar in that their reaction rates are faster at higher substrate concentrations, slower at higher product con-centrations, and not affected by other metabolites.
Two-way monotonic kinetics [79, 84] is defined as:
Bvjpxq Bxi $ ' & ' %
¡ 0, if xi is a substrate concentration of reaction j
¤ 0, if xi is a product concentration of reaction j
0, otherwise.
(2.11)
CRNs with this property are referred to as having ’two-way mono-tonic kinetics’. The literature also refers to this type of kinetics as
tonic kinetics is very general. For example, mass action, Michaelis-Menten, and Hill are two-way monotonic.
In systems with two-way monotonic kinetics there is an intimate relationship between the stoichiometric component, represented by sto-ichiometric matrix S and the kinetic component vpxq: the sign pattern of Sij P Z determines the sign of BvBxjpxq
i .
Note that the assumption of two-way monotonic kinetics excludes the possibility of regulation. For example, it is not possible for strates to inhibit the reaction or for metabolites other than the sub-strates or products to influence the reaction rate.
2.3
Dynamical properties of CRNs
Most results of CRNT concern characterization of steady states and their stability. Other dynamical properties have also been investigated extensively, such as persistence, permanence, and boundedness, but this thesis will only address the results concerning steady-state prop-erties.
Equilibrium represents the state of a system when all concentrations of reactants remain constant in time. Mathematically speaking, an equilibrium x is a solution to the equation (2.4) set to zero (9x 0). Alternatively,
Svpxq 0. (2.12)
The concepts of equilibrium and steady state are not identical, although they are closely related. If a chemical reaction is closed, x is a vector of equilibrium concentrations, or simply an equilibrium. If a CRN is open, x is called a steady state.
A steady state does not always exist, but when it does, it can exist as a unique steady state or multiple steady states (multistationarity). A steady state is locally asymptotically stable if there exists a local neighborhood such that if the initial condition is in this neighborhood, then the system will eventually converge to this steady state [87].
A steady state is globally asymptotically stable if all solutions ap-proach the steady state from any point of state space [87].
The steady-state properties and their stability are investigated with the use of the Jacobian matrix. The Jacobian matrix is defined as:
Jpxq B 9x Bx BSvpxq Bx S Bvpxq Bx . (2.13)
If all eigenvalues of the Jacobian matrix of a given system have negative real parts, then the steady state in such a system is locally asymptotically stable [88].
2.4
General kinetics
The structure of this section is based on the theoretical approaches that identify the number of steady states and their stability. The first approach is based on the properties of stoichiometric and Jaco-bian matrices. The second approach is based on the properties of the species-reaction graph.
2.4.1
Stoichiometric and Jacobian matrices
As I discussed when describing two-way monotonic kinetics, the sign pattern of Sij determines the sign of BvBxjpxq
i . As a result, the sign pattern of the stoichiometric matrix determines the sign pattern of the Jaco-bian matrix (equation 2.13). This is a very powerful insight, because in many cases it is enough to know the sign pattern of the Jacobian matrix to prove its general properties, in particular the properties of eigenvalues.
Using sign patterns, Banaji and Beigent [77] proved that in net-works with single-substrate-single-product reactions with binary stoi-chiometry and two-way monotonic kinetics, the Jacobian matrix is a negative M-matrix, i.e., a matrix with non-negative off-diagonal en-tries and all leading principal minors alternating in sign. Negative M-matrices have eigenvalues with negative real parts [89]. Therefore, networks for which Jacobian matrices are negative M-matrices admit unique and globally asymptotically stable steady states.
Using a similar approach, in Chapter 3 I prove uniqueness and local asymptotic stability for networks with single-substrate-single-product
vided that there are no cycles; for networks with metabolic cycles, I show the condition for local asymptotic stability of steady states [79]. In Chapter 3 I also show that if, in metabolic networks with multiple-substrate-multiple-product reactions and with monotonic kinetics, the stoichiometric matrix is singular, then steady states are in the form of manifolds and the dimension of manifold of steady states is larger than or equal to the rank of the stoichiometric matrix. Moreover, I conjec-ture that these manifolds of steady states are locally asymptotically stable [79].
Banaji et al. [86] defined a ’sign-nonsingular’ matrix as a matrix with a nonzero determinant that can be determined from the signs of the matrix entries. A matrix is ’strongly sign determined’ (SSD) if its submatrices are either singular or sign-nonsingular [86]. Therefore, CFSTR networks with two-way monotonic kinetics and SSD properties cannot admit multiple steady states [86]. That is, if the steady state exists, it is unique.
2.4.2
Species-reaction graph
A species-reaction graph (SR-graph) is a bipartite undirected graph with chemical species and reactions as node types. A chemical species node is connected to the reaction node if it participates in the cor-responding reaction. Every edge in the SR-graph has a weight that corresponds to the stoichiometric coefficient with which the chemical species participates in a certain reaction. Moreover, every edge has a symbol of a reaction complex to which a corresponding species belongs (Figure 2.1).
A complex pair (c-pair) is a pair of edges in the SR-graph that are connected to the same reaction node and carry the same complex label. For example, in the SR-graph in Figure 2.1 that corresponds to the reaction network (2.1), the c-pairs (color marked) with the cor-responding complex labels are: red (E S), green (F P ), orange (E P ), and blue (E S).
A cycle is a closed path in which no node or edge is traversed twice. There are six cycles in the SR-graph in Figure 2.1. An odd-cycle
(o-Figure 2.1: Species-reaction graph for network (2.1). Circles represent chemical species nodes. Rectangles represent reaction nodes. Colors of edges represent different closed pairs (c-pairs) and labels next to edges represent different complexes.
cycle) in the SR-graph contains an odd number of c-pairs, whereas an even-cycle (e-cycle) contains an even number of c-pairs. In the SR-graph in Figure 2.1, the largest outer cycle contains all four c-pairs; it is therefore an e-cycle. A stoichiometric cycle is one for which alternately multiplying and dividing the stoichiometric coefficients along the cycle gives the final result of one. A 1-cycle is a cycle in which every weight associated with every edge is 1, assuming all stoichiometric coefficients are 1, which applies to most biochemical networks [76]. Note that a 1-cycle is a special case of an s-cycle. Two cycles have a species-to-reaction (S-to-R) intersection if their common edges constitute a path that begins at a species node and ends at a reaction node or if they constitute a disjoint union of such paths. If the path of common edges begins at a reaction node and ends at a reaction node, the two cycles have a reaction-to-reaction (R-to-R) intersection.
Now the main results of CRNT with respect to SR-graphs can be formulated.
Theorem 1 Consider a CFSTR CRN such that in its SR-graph: • each cycle is an o-cycle, an s-cycle, or both
• no two e-cycles have an S-to-R intersection
then, the reaction network does not have the capacity for multiple pos-itive steady states for mass action [90, 91], weakly monotonic [92], or two-way monotonic [85] kinetics.
The requirement of the fully open network (CFSTR) is dropped by the introduction of the concept of ’entrapped species’, i.e., reactants that do not flow into or out of the network. In this case, it was shown that Theorem 1 holds in the case of weakly reversible networks with mass action kinetics [83] and those with arbitrary kinetics for non-degenerate steady states (i.e., for which the Jacobian matrix at the steady state x is nonzero, Jpxq 0)[82].
As an application of this theorem, consider the SR-graph in Figure 2.1. In this case, all cycles are e-cycles. The weights of all edges are equal to one; thus, all cycles are also 1-cycles. Moreover, all intersec-tions of e-cycles are of the R-to-R type. Therefore, the condiintersec-tions of Theorem 1 are fulfilled and reaction network (2.1), assuming it is a CFSTR with, for example, two-way monotonic kinetics, does not have the capacity for multiple steady states, irrespective of its parameters.
In CRNs with directed, strongly connected SR-graphs (i.e., graphs for which there are directed paths from every node to every other node) and NAC kinetics combined with the special partitioning of the stoi-chiometric matrix, the equilibrium is unique and locally asymptotically stable [93]. If, in addition, all reactions are reversible and the steady state is positive, then it is also globally asymptotically stable [93].
2.5
Mass action kinetics
More precise results can be derived if one is willing to make specific assumptions regarding the kinetics of a system. In fact, most of the research on CRNT has focused on a mathematical analysis of systems with mass action kinetics. Mass action kinetics is a specific case of
the general kinetics types discussed above. Accordingly, all the results discussed thus far apply to mass action systems as well.
As with the previous section, the structure of this section is based on theoretical approaches that identify the number of CRN steady states and their stability. The first approach is based on the properties of the stoichiometric and Jacobian matrices. The second approach is based on the concept of deficiency, and the third approach is based on concepts of balancing.
Before presenting the main results of models with mass action ki-netics, I will first introduce the concepts of stoichiometric subspace and stoichiometric compatibility class.
A linear combination of columns of the stoichiometric matrix, S spanpSq, is called the stoichiometric subspace [94]. The number of linearly independent columns of the stoichiometric matrix, s, is equal to the dimension of the stoichiometric subspace, s dimpSq. Alter-natively, the dimension of the stoichiometric subspace is equal to the rank of the stoichiometric matrix. The dynamics of chemical species concentrations occur in the stoichiometric subspace.
For example, there are six metabolites in the futile cycle (2.1). However, the number of linearly independent columns of the stoichio-metric matrix (equation 2.3) is three. Therefore, the dynamics is con-strained to the three-dimensional subspace, rather than the possible six-dimensional space. The other three dimensions correspond to con-servation laws for this system. Specifically, F F P const, E ES const, and S ES P F P const are conserved quantities.
A stoichiometric compatibility class (SCC) is a translation of the stoichiometric subspace determined by the initial condition SCpx0q
x0 S. The solutions of equations that describe the concentration
dynamics of the reactants in CRNs cannot occur in the whole space of possible concentrations, but are restricted to stoichiometric compati-bility classes. That is, given the initial condition x0, the concentration
dynamics of the reactant occur only in the corresponding SCC, SCpx0q.
For example, consider the following CRN
0 5 10 x1 20 30 40 50 x2 0 10 20 x3
Figure 2.2: Stoichiometric compatibility classes for the network (2.14). Parameters: k1 1, k2 0.1. Initial conditions x1p0q 10, x2p0q 25 (blue lines), x2p0q 50 (red lines), x3p0q is from 5 to 15 with a step size of one.
In this network, the quantity x1 x2 const is conserved.
There-fore, the rank of the stoichiometric matrix is equal to two. The number of linearly independent columns of the stoichiometric matrix is also equal to two. Therefore, the dimension of the stoichiometric subspace is equal to two as well. Hence, the metabolite concentration dynamics do not occur in three-dimensional space (as do the number of chemi-cal species), but occur in a two-dimensional plane. To illustrate this, I assume the following set of differential equations that describe the dynamical behavior of network (2.14):
9x1 9x2 k1x1x2
9x3 k1x1x2 k2x3.
(2.15) The solutions for the system of equations (2.15) are presented in Figure 2.2. Blue and red lines belong to different stoichiometric com-patibility classes (SCCs) in which the dynamics of metabolite concen-trations occur. These SCCs differ by the initial condition x2p0q. The
blue SCC corresponds to x2p0q 25 and the red SCC corresponds to
x2p0q 50.
2.5.1
Stoichiometric and Jacobian matrices
For CFSTR networks with mass action kinetics, Craciun and Feinberg [95] showed that a network has the capacity for multiple steady states only if the determinant of the Jacobian matrix equals zero, but not if it does not equal zero. Banaji and Craciun [85] showed that for CFSTR networks with mass action kinetics, the SSD property of the stoichio-metric matrix and properties of SR-graph that ensure impossibility of multiple steady states (Theorem 1) are equivalent.
2.5.2
Balancing
The concepts of thermodynamic balancing (detailed balancing) and complex balancing are key to many branches of CRNT.
An equilibrium x is a thermodynamic equilibrium if vpxq 0, i.e., the forward reaction rate is equal to the reverse reaction rate for every reaction. Any thermodynamic equilibrium is an equilibrium (in the sense of equation 2.12), but not all equilibriums are thermodynamic in nature. A CRN that admits thermodynamic equilibrium is thermody-namically balanced. For example, consider the reaction X1
v
é X2 with
vpxq k1x1 k2x2. If this reaction is thermodynamically balanced,
then k1x1 k2x2. In the literature, thermodynamic balancing is also
referred as detailed balancing [78].
Van der Schaft et al. proved that closed and thermodynamically balanced CRNs with mass action kinetics have global asymptotic con-vergence to a set of equilibria (typically not a unique equilibrium) [96]. If a closed network with mass action kinetics is thermodynam-ically balanced [97], then there is a unique equilibrium in each SCC [98, 99, 100, 101]. If the CRN is closed, the total number of SCCs may be infinite as may the number of equilibria. However, there is a unique equilibrium for each SCC (i.e., each initial condition), which precludes complex types of dynamical behavior, e.g., multistability or oscillations.
namic balancing [99]. In a complex balanced system, the net flow to each reaction complex equals the net flow out of the complex [99]. Chemically, this means that for the complex balanced equilibrium x, the concentrations of all chemical species and complexes remain con-stant [78]. Complex balanced networks with mass action kinetics have several strong properties [98, 99, 100, 101]:
1. Every complex balanced network is weakly reversible.
2. In a mass action system that permits a complex balanced equi-librium, all equilibria must be complex balanced.
3. There is a unique equilibrium in each SCC, i.e., for each initial condition.
4. The equilibrium in each SCC is locally asymptotically stable. Therefore, in complex balanced networks with mass action kinetics, dynamical behavior such as multistability, oscillations, and chaos are not possible [99].
In 1974, Horn stated the conjecture [102], which was later named the Global Attractor Conjecture [103]. The Global Attractor Conjec-ture was proved only recently [104], forty years after its introduc-tion, although results for special cases were obtained prior to this [103, 105, 106, 107, 108, 109, 110].
Global Attractor Theorem: the complex balanced equilib-rium in a system with mass action kinetics is globally asymptotically stable for its SCC.
Another result states that, in an open complex-balanced CRNT with constant inflows and proportional outflows with mass action kinet-ics, if the complex stoichiometric matrix is surjective, then the steady state is unique and locally asymptotically stable, given mild additional assumptions [78].
2.5.3
Deficiency
δ m l s. (2.16) Here, m is the number of different reaction complexes, l is the number of different linkage classes, and s is the dimension of the stoichiometric subspace, which is equal to the number of linearly independent columns of the stoichiometric matrix.
The first result of CRNT was the ’deficiency zero theorem’ [94, 99, 98, 100]. In the case of zero-deficiency, there is a one-to-one correspon-dence between the rate vector 9x of chemical species x and the rate vector 9y of complexes y [96].
Deficiency Zero Theorem When considering a mass action network of zero deficiency, the following statements hold true for any parameters:
1. if the network is not weakly reversible, then for arbitrary kinetics (e.g., mass action, monotonic, etc.), the system cannot have a positive equilibrium or sustained oscillations.
2. if the network is weakly reversible, then for mass action kinetics, the network has a single positive equilibrium point in each SCC, and that equilibrium point is locally asymptotically stable.
Consider the following CRN (Example 3.4.2 in [111]): A1 ÝÑ A2
Ô Ö 2A1 ÐÑ 2A3.
A3 A4
(2.17)
This system has five reaction complexes: tA1, A2, A3 A4, 2A1, 2A3u,
two linkage classes, l 2, and it is weakly reversible. The stoichiomet-ric matrix for this system is:
S 1 0 1 2 1 1 0 0 0 1 1 2 0 1 1 0 . (2.18)
The number of linearly independent columns in this matrix is 3. The deficiency of the network is therefore δ n l s 5 2 3 0.
steady equilibrium for each SCC.
A strongly connected network component is a subnetwork in which there is a direct path from each node to every other node. A strongly connected network component is called terminal if there is no reaction leading out of the set. Consider the following example:
C1 Ø C2 Ñ C3 Ø C4. (2.19)
There are two strongly connected sets of reaction complexes C11 tC1, C2u and C21 tC3, C4u. However, only the complex set C21 is
terminal.
Now we can formulate the second classic result of CRNT, the ’de-ficiency one theorem’, which was proved by M. Feinberg [94].
Deficiency One Theorem: Consider a CRN with mass action kinetics with deficiency δ and l linkage classes. Let δi, i 1, ..., l denote
the deficiencies of the individual linkage classes considered as their own networks. Suppose the following is true:
1. δi ¤ 1 for all i 1, ..., l,
2. °li1δi δ
3. each linkage class contains a single terminal strongly-linked com-ponent.
Then, if such a system admits a positive equilibrium, there exists precisely one equilibrium in each positive SCC. Furthermore, if that network is weakly reversible, every mass action system permitted by the network has a positive equilibrium.
Network (2.1) has n 6 reaction complexes, l 2 linkage classes, and s 3 dimensions of the stoichiometric subspace. Hence, the net-work deficiency is δ 1. The deficiencies of the linkage classes are: δ1 δ2 δ. Therefore, the deficiency one theorem may not be applied
2.6
Concluding remarks
CRNT is a broad area that can be developed in many different di-rections; however, various parallel directions have recently begun to merge. An overview of the results of CRNT discussed in this chapter is provided in Table 1.
CRNT has natural applications in biology. For example, knowledge gained and techniques developed from CRNT were applied to analyze metabolic networks [79, 80, 81, 109], mitochondrial electron transport [77], T-cell receptor signal transduction [110, 112], multisite protein phosphorylation [75, 113, 114], to model anti-cancer drugs effects [91] and interactions of interleukin-1 receptors with antagonists [115].
Interestingly, different approaches of CRNT predict that in CRNs steady states are unique and locally or globally asymptotically stable under general assumptions on network structure and kinetics, irrespec-tively of parameter values [79, 116]. The result is that, in such net-works, dynamical behavior is very simple: monotonic convergence to a unique steady state or to a set of steady states. Complex dynamical behavior, such as oscillations, multistability, or chaos, are excluded.
CRNT is a powerful tool for studying qualitative behavior of bio-chemical networks in the presence of structural information only. The generality of CRNT introduces its own disadvantages. Specifically, CRNT is not directly applicable if one is interested in specific trajec-tories or specific values of steady states. In this case, detailed modeling, often combined with numerical simulations, should be used. Neverthe-less, CRNT can provide valuable support in such situations by pre-dicting the type of dynamical behavior to expect (and not to expect) for all choices of parameters.
Structure Balancing Deficiency In-/outflo w Mixed Results Remarks Refs Kinetics SR SSD SSSP MSMP BS WR TB CB D0 D1 Cl CFSTR ES I/O J U M L G Arb + -+ [79] + + + 1 [82] NA C + + + + + [77] + + + + 2 [79] + + + [86] + + + [85] + + 3 [93] + + 4 [93] WM + + + [92] MA + + + [96] + + + + scc [97] + + + + scc [98] + + + + 5 [78] + + + + + scc [99] + + + + + scc, 6 [94] + + + + [83] + + + [90] + + + [95] + -+ [95] T able 2.1 : Summary of CRNT results discussed in this chapter. Assumptions abbreviations : Arb -arbitrary kinetics, NA C -non-auto catalytic (monotonic, N1C) kinetics; WM -w e a kly monotonic kinetics; SR -con di tions on SR-gr aph from Theorem 1; SSD -strongly sign determined sto ic hiometric matrix; SSSP -single-substrate-single-pro duct net w ork; MSMP -m ultiple-substrate-m ultiple-pro duct net w ork; BS -binar y stoic hiometry; WR -w eak rev ersibilit y; T B -thermo dynamic balancing; CB -complex balancing; D0 -zero deficiency; D1 -conditi o ns of deficiency from deficiency one theorem; Cl -closed n e tw ork; CF STR -con tin uous flo w stirred-tank reactor; ES -e n trapp ed sp ecies; I/O -net w o rk has inflo ws and outflo ws of some reactan J -Jacobian matrix is nonsingular. Results abbreviations : U -unique steady state or equilibrium (if one exists); M -m ul tiple steady states; L -lo asymptotic sta bi lit y o f steady state; G -global asymptotic stabilit y of steady state; Remarks : scc -for a stoic hiometric compatibilit y class (i.e., for giv en initial condition); 1 -steady states are non-degenerativ e (i .e., if the determinan t of the Jacobian matrix is nonzero); 2 -the net w or k has no cycles; a directed sp ecies-reaction graph is strongly connected with sp ecial partitioning of the stoic hiometric matrix; 4 -in addition to a ss u m ption 3, steady state is p ositiv e and the net w ork is rev ersible; 5 -compl ex stoic hiometric matrix is surjectiv e and lo cal asymptotic sta bi lit y is pro vided under the assumption of p ersistence (if ev ery chemical sp ecies is presen t at the initial time, none of the chemical sp ecies reac hes zero concen tration in the course of time); eac h link age class con ta ins a single terminal strongly-link ed comp onen t.
Steady states and their stability in
metabolic networks without regulation
Abstract
Metabolic networks are often extremely complex. Despite intensive efforts many details of these networks, e.g., ex-act kinetic rates and parameters of metabolic reex-actions, are not known, making it difficult to derive their prop-erties. Considerable effort has been made to develop a theory about properties of steady states in metabolic net-works that are valid for any values of parameters. General results on the uniqueness of steady states and their sta-bility have been derived with specific assumptions on re-action kinetics, stoichiometry and network topology. For example, deep results have been obtained under the as-sumptions of mass action reaction kinetics, continuous flow stirred tank reactors (CFSTR), concordant reaction net-works and others. Nevertheless, a general theory about properties of steady states in metabolic networks is still missing. Here we make a step further in the quest for such a theory. Specifically, we study properties of steady states in metabolic networks with monotonic kinetics in re-lation to their stoichiometry (simple and general) and the number of metabolites participating in every reaction (sin-gle or many). Our approach is based on the investigation of properties of the Jacobian matrix. We show that
sto-ichiometry, network topology and the number of metabo-lites that participate in every reaction have a large influ-ence on the number of steady states and their stability in metabolic networks. Specifically, metabolic networks with single-substrate-single-product reactions have disconnected steady states, whereas in metabolic networks with multiple-substrates-multiple-product reactions manifolds of steady states arise. Metabolic networks with simple stoichiometry have either a unique globally asymptotically stable steady state or asymptotically stable manifolds of steady states. In metabolic networks with general stoichiometry the steady states are not always stable and we provide conditions for their stability. To demonstrate the biological relevance we illustrate the results on the examples of the TCA cycle, the mevalonate pathway and the Calvin cycle.
