1. Computational and Systems Biology
Download icon

Emergence and propagation of epistasis in metabolic networks

  1. Sergey Kryazhimskiy  Is a corresponding author
  1. Division of Biological Sciences, University of California, San Diego, United States
Research Article
  • Cited 0
  • Views 1,509
  • Annotations
Cite this article as: eLife 2021;10:e60200 doi: 10.7554/eLife.60200
Voice your concerns about research culture and research communication: Have your say in our 7th annual survey.

Abstract

Epistasis is often used to probe functional relationships between genes, and it plays an important role in evolution. However, we lack theory to understand how functional relationships at the molecular level translate into epistasis at the level of whole-organism phenotypes, such as fitness. Here, I derive two rules for how epistasis between mutations with small effects propagates from lower- to higher-level phenotypes in a hierarchical metabolic network with first-order kinetics and how such epistasis depends on topology. Most importantly, weak epistasis at a lower level may be distorted as it propagates to higher levels. Computational analyses show that epistasis in more realistic models likely follows similar, albeit more complex, patterns. These results suggest that pairwise inter-gene epistasis should be common, and it should generically depend on the genetic background and environment. Furthermore, the epistasis coefficients measured for high-level phenotypes may not be sufficient to fully infer the underlying functional relationships.

Introduction

Life emerges from an orchestrated performance of complex regulatory and metabolic networks within cells. The blueprint for these networks is encoded in the genome. Mutations alter the genome. Some of them, once decoded by the cell, perturb cellular networks and thereby change the phenotypes important for life. Understanding how mutations affect the function of cellular networks is key to solving many practical and fundamental problems, such as finding mechanistic causes of genetic disorders (Hu et al., 2011; Fang et al., 2019), deciphering the architecture of complex traits (Zuk et al., 2012; Mackay, 2014; Wei et al., 2014), building artificial cells (Hutchison et al., 2016), explaining past, and predicting future evolution (Blount et al., 2008; Wiser et al., 2013; de Visser and Krug, 2014; Harms and Thornton, 2014; Kryazhimskiy et al., 2014; Sailer and Harms, 2017a; Sohail et al., 2017). Conversely, mutations can help us learn how cellular networks are organized (Phillips, 2008; van Opijnen and Camilli, 2013).

To infer the wiring diagram of a cellular network that produces a certain phenotype, one approach in genetics is to measure the pairwise and higher-order genetic interactions (or ‘epistasis’) between mutations that perturb it (Phillips, 2008). Much effort has been devoted in the past 20 years to a systematic collection of such genetic interaction data for several model organisms and cell lines (Kelley and Ideker, 2005; Lehner et al., 2006; Jasnos and Korona, 2007; Collins et al., 2007; St Onge et al., 2007; Typas et al., 2008; Roguev et al., 2008; Costanzo et al., 2010; Szappanos et al., 2011; Huang et al., 2012; Roguev et al., 2013; Bassik et al., 2013; Babu et al., 2014; Costanzo et al., 2016; van Leeuwen et al., 2016; Skwark et al., 2017; Du et al., 2017; Heigwer et al., 2018; Horlbeck et al., 2018; Norman et al., 2019; Liu et al., 2019; Kuzmin et al., 2018; New and Lehner, 2019; Celaj et al., 2020). This approach is particulary powerful when the phenotypic effect of one mutation changes qualitatively depending on the presence or absence of a second mutation in another gene, for example when a mutation has no effect on the phenotype in the wildtype background, but abolishes the phenotype when introduced together with another mutation, such as synthetic lethality (Tong et al., 2001). Such qualitative genetic interactions can often be directly interpreted in terms of a functional relationship between gene products (Tong et al., 2001; Davierwala et al., 2005; Phillips, 2008).

Most pairs of mutations do not exhibit qualitative genetic interactions. Instead, the phenotypic effect of a mutation may change measurably but not qualitatively depending on the presence or absence of other mutations in the genome (Babu et al., 2014; Costanzo et al., 2016). The genetic interactions can in this case be quantified with one of several metrics that are termed ‘epistasis coefficients’ (Wagner et al., 1998; Hansen and Wagner, 2001; Mani et al., 2008; Wagner, 2015; Poelwijk et al., 2016). Although some rules have been proposed for interpreting epistasis coefficients, in particular, their sign (Dixon et al., 2009; Lehner, 2011; Baryshnikova et al., 2013), the validity, and robustness of these rules are unknown because there is no theory for how functional relationships translate into measurable epistasis coefficients in any system (Lehner, 2011; Domingo et al., 2019). To avoid this major difficulty, most large-scale empirical studies focus on correlations between epistasis coefficients rather than on their actual values (but see Velenich and Gore, 2013, for a notable exception). Genes with highly correlated epistasis profiles are then interpreted as being functionally related (Segrè et al., 2005; Bellay et al., 2011; Babu et al., 2014; Costanzo et al., 2016; Horlbeck et al., 2018). Although this approach successfully groups genes into protein complexes and larger functional modules (Michaut et al., 2011; Bellay et al., 2011), it does not reveal the functional relationships themselves. As a result, many if not most, genetic interactions between genes and modules still await their biological interpretation (Costanzo et al., 2016; Fang et al., 2019).

While geneticists measure epistasis to learn the architecture of biological networks, evolutionary biologists face the reverse problem: they need to know how the genetic architecture constrains epistasis at the level of fitness. Epistasis determines the structure of fitness landscapes on which populations evolve (Fragata et al., 2019). Understanding it would bear on many important evolutionary questions, such as why so many organisms reproduce sexually (Kondrashov, 2018), how novel phenotypes evolve (Blount et al., 2008; Bridgham et al., 2009; Natarajan et al., 2013; Harms and Thornton, 2014), how predictable evolution is (Weinreich et al., 2006; Tenaillon et al., 2012; Wiser et al., 2013; Kryazhimskiy et al., 2014), etc. So far, evolutionary biologists have relied primarily on abstract models of fitness landscapes (see Orr, 2005, for a review), rather than those firmly grounded in organismal biochemistry and physiology (e.g. Dykhuizen et al., 1987; Das et al., 2020). For example, Fisher’s geometric model—one of the most widely used fitness landscape models—is explicitly devoid of the physiological and biochemical details (Fisher, 1930; Tenaillon, 2014; Martin, 2014).

A theory of epistasis must address two challenges. First, it must specify how the architecture of a biological network constrains epistasis. Such knowledge is important not only for evolutionary questions, but also for the inference problem in genetics. Consider a biological network module that produces a phenotype of interest but whose internal structure is unknown. By genetically perturbing all genes within the module and measuring the phenotype in all single, double and possibly some higher-order mutants, we can obtain the matrix of epistasis coefficients. In principle, we can then fit a network topology and parameters to these data. However, without knowing what information about the network is contained in the matrix in the first place, we cannot be sure whether the inferred topology and parameters are close to their true values or represent one of many possible solutions consistent with the data.

The second challenge is that epistasis may arise at a different level of biological organization than where it is measured by the experimentalist or by natural selection. For example, geneticists are often interested in understanding the structures of specific regulatory or metabolic network modules (Collins et al., 2007; Costanzo et al., 2010). However, measuring the peformance of a module directly is often experimentally difficult or impossible. Then epistasis is measured for an experimentally accessible ‘high-level’ phenotype, such as fitness, which depends on the performance of the focal ‘lower-level’ module, but also on other unrelated modules. However, if we do not know how epistasis that originally emerged in one module maps onto epistasis that is measured, it is unclear what we can infer about module’s internal structure.

Evolutionary biologists encounter a related problem when they wish to learn the evolutionary history of a protein or a larger cellular module. To do so, they would in principle need to know how different mutations in this module affected fitness of the whole organism in its past environment. But such information is rarely available. Instead, it is sometimes possible to reconstruct past mutations and measure their biochemical effects in the lab (Lunzer et al., 2005; Bridgham et al., 2009; Natarajan et al., 2013; Sarkisyan et al., 2016). When interesting patterns of epistasis are identified at the biochemical level, it is usually assumed that the same patterns manifested themselved at the level of fitness and drove module’s evolution. However, this is not obvious. If interactions with other modules distort epistasis as it propagates from the biochemical level to the level of fitness (Snitkin and Segrè, 2011), our ability to infer past evolutionary history from in vitro biochemical measurements could be diminished. Therefore, the second challenge that a theory of epistasis must address is how epistasis propagates from lower-level phenotypes to higher-level phenotypes.

There is a large body of theoretical and computational literature on epistasis. As early as 1934, Sewall Wright realized that epistasis naturally emerges in molecular networks (Wright, 1934). This was later explicitly demonstrated in many mathematical and computational models (e.g. Kacser and Burns, 1981; Keightley, 1989; Szathmáry, 1993; Gibson, 1996; Keightley, 1996; Wagner et al., 1998; Omholt et al., 2000; Peccoud et al., 2004; Gjuvsland et al., 2007; Gertz et al., 2010; Fiévet et al., 2010; Pumir and Shraiman, 2011). Metabolic control analysis became one of the most successful and general frameworks for understanding epistasis between metabolic genes (Kacser and Burns, 1973). Dean et al., 1986, Dykhuizen et al., 1987, Dean, 1989, Lunzer et al., 2005, MacLean, 2010 used it to interpret the empirically measured fitness effects of mutations and their interactions in terms of the metabolic relationships between the products of mutated genes. Kacser and Burns, 1981, Hartl et al., 1985, Keightley, 1989, Clark, 1991, Keightley, 1996, Bagheri-Chaichian et al., 2003, Bagheri and Wagner, 2004, Fiévet et al., 2010 explored the implications of epistasis in metabolism for genetic variation in populations, their response to selection, long-term evolutionary dynamics and outcomes, such as the evolution of dominance. However, most studies analyzed only the linear metabolic pathway (but see Keightley, 1989) and assumed that fitness equals flux through the pathway (but see Szathmáry, 1993), thereby bypassing the problem of epistasis propagation.

There have been few attempts to theoretically relate the molecular architecture of an organism to the types of epistasis that would arise for its high-level phenotypes, such as fitness. Segrè et al., 2005 and He et al., 2010 used flux balance analysis (FBA, Orth et al., 2010) to compute genome-wide distributions of epistasis coefficients in metabolic models of Escherichia coli and Saccharomyces cerevisiae and arrived at starkly discordant conclusions. Recently, Alzoubi et al., 2019 showed that FBA is generally poor in predicting experimentally measured genetic interactions, suggesting that it might be difficult to understand the emergence and propagation of epistasis by relying exclusively on genome-scale computational models. Sanjuán and Nebot, 2008 and Macía et al., 2012 modeled various abstract metabolic and regulatory networks and found a possible link between epistasis and network complexity. The work by Chiu et al., 2012 is a more systematic attempt to develop a general theory of epistasis. They established a fundamental connection between epistasis and the curvature of the function that maps lower-level phenotypes onto a higher-level phenotype. However, further progress has been so far hindered by uncertainty in what types of functions map phenotypes onto one another in real biological systems. Previous studies made various idiosyncratic choices with respect to this mapping, leaving us without a clear guidance as to the conditions or systems where they are expected to hold.

To overcome this problem, here I consider a whole class of hierarchical metabolic networks and obtain the family of all functions that determine how the effective activity of a larger metabolic module can depend on the activities of smaller constituent modules. There are several advantages to this approach. First, it leads to an intuitive understanding of how the structure of the network influences epistasis emergence and propagation. Second, my approach is based on basic biochemical principles, so it should be relevant for many phenotypes. For example, epistasis is often measured at the level of growth rate (Jasnos and Korona, 2007; St Onge et al., 2007; Babu et al., 2014; Costanzo et al., 2016), and metabolism fuels growth. Moreover, metabolic genes occupy a large fraction of most genomes (Orth et al., 2011) and the general organization of metabolism is conserved throughout life (Csete and Doyle, 2004). Thus, by understanding genetic interactions between metabolic genes, we will gain an understanding of a large fraction of all genetic interactions.

In my model, I consider a hierarchical network with first-order kinetics but arbitrary topology, and ask two questions related to the two challenges mentioned above. (1) How does an epistasis coefficient that arose at some level of the metabolic hierarchy propagate to higher levels of the hierarchy? (2) How does the network topology constrain the value of an epistasis coefficient between two mutations that affect different enzymes in this network? I obtain answers to these questions analytically in the limiting case when the effects of mutations are vanishingly small. I then computationally probe the validity of the conclusions outside of the domain where they are expected to hold.

My model is not intended to generate predictions of epistasis for any specific organism. Instead, its main purpose is to provide a baseline expectation for how epistasis that emerges at lower-level phenotypes manifests itself at higher-level whole-organism phenotypes, such as fitness, and what kind of information may be gained from measurements of such higher-level epistasis. One possible outcome of this analysis is that there may be fundamental limitations to what an epistasis measurement at one level of biological organization can tell us about epistasis at another level. On the other hand, if it turns out that there is a general correspondence between epistasis coefficients at different levels in this simple model, then it may be worth developing more sophisticated and general models on which inference from data can be based.

Model

Hierarchical metabolic network

Consider a set of metabolites A={1,2,,n} with concentrations S1,,Sn which can be interconverted by reversible first-order biochemical reactions. The rate of the reaction converting metabolite i into metabolite j is xij(Si-Sj/Kij) where Kij is the equilibrium constant. The rate constants xij, which satisfy the Haldane relationships xji=xij/Kij (Cornish-Bowden, 2013), form the matrix x=xiji,j=1n. The metabolite set A and the rate matrix x define a biochemical network 𝒩=(A,x).

The first-order kinetics assumption makes the model analytically tractable, as discussed below; biologically, it is equivalent to assuming that all enzymes are far from saturation. The rate constants xij depend on the concentrations and the specific activities of enzymes and therefore can be altered by mutations. Kij characterize the fundamental chemical nature of metabolites i and j and cannot be altered by mutations (Savageau, 1976).

The whole-cell metabolic network is large, and it is often useful to divide it into subnetworks that carry out certain functions important for the cell. I define subnetworks mathematically as follows. I say that two metabolites i and j are adjacent (in the graph-theoretic sense) if there exists an enzyme that catalyzes a biochemical reaction between them, that is, if xij>0. Now consider a subset of metabolites AμA. For this subset, let AμIO be the set of all metabolites that do not belong to Aμ but are adjacent to at least one metabolite from Aμ. Let xμ be the submatrix of x which corresponds to all reactions where both the product and the substrate belong to either Aμ or AμIO. The metabolite subset Aμ and the rate matrix xμ form a subnetwork μ=(Aμ,xμ) of network 𝒩. I refer to Aμ and AμIO as the sets of internal and ‘input/output’ (‘I/O’ for short) metabolites for subnetwork μ, respectively. Thus, all internal metabolites and all reactions that involve only internal and I/O metabolites are part of the subnetwork. Note that the I/O metabolites do not themselves belong to the subnetwork, but reactions between them, if they exist, are part of the subnetwork. Metabolites that are neither internal nor I/O for μ are referred to as external to subnetwork μ. These definitions are illustrated in Figure 1A.

Illustration of a hierarchical metabolic network and its coarse-graining.

(A) White rectangle represents the whole metabolic network 𝒩. Example subnetworks μ and ν are represented by the dark and light gray rectangles. Only metabolites and reactions that belong to these subnetworks are shown; other metabolites and reactions in 𝒩 are not shown. Metabolites 1 and 5 may be adjacent to other metabolites in 𝒩; this fact is represented by short black lines that do not terminate in metabolites. Subnetworks μ and ν are both modules because there exists a simple path connecting their I/O metabolites that lies within μ and ν and contains all their internal metabolites (dashed blue line). (B) Network 𝒩 can be coarse-grained by replacing module μ at steady state with an effective reaction between its I/O metabolites 1 and 2, with the rate constant is yμ. (C) Network 𝒩 can be coarse-grained by replacing module ν at steady state with an effective reaction between its I/O metabolites 1 and 5, with the rate constant is yν.

The main objects in this work are biochemical modules, which are a special type of subnetworks. To define modules, I introduce some auxiliary concepts. I say that two metabolites i and j are connected if there exists a series of enzymes that interconvert i and j, possibly through a series of intermediates. Mathematically, i and j are connected if there exists a simple (i.e. non-self-intersecting) path between them. If all metabolites in this path are internal to the subnetwork μ (possibly excluding the terminal metabolites i and j themselves) then i and j are connected within the subnetwork μ, and such path is said to lie within μ. By this definition, metabolites i and j can be connected within μ only if they are either internal or I/O metabolites for μ.

Definition 1

A subnetwork μ is called a module if (a) it has two I/O metabolites, and (b) for every internal metabolite iAμ, there exists a simple path between the I/O metabolites that lies within μ and contains i.

This definition is illustrated in Figure 1A. The assumption that modules only have two I/O metabolites is not essential. However, mathematical calculations become unwieldy when the number of I/O metabolites increases. Moreover, modules with just two I/O metabolites already capture two most salient features of metabolism: its directionality, and its complex branched topology (Csete and Doyle, 2004). Such modules are a natural generalization of the linear metabolic pathway which has been extensively studied in the previous literature (Kacser and Burns, 1973; Szathmáry, 1993; Bagheri-Chaichian et al., 2003; MacLean, 2010).

Modules have two important properties. First, for any given concentrations of the two I/O metabolites, all internal metabolites in the module can achieve a unique steady state which depends only on concentrations of these I/O metabolites but not on the concentrations of any other metabolites in the network (see Proposition 1 in Materials and methods). Now consider a module μ whose I/O metabolites are (without loss of generality) labeled 1 and 2 (Figure 1A). The second property is that, at steady state, the flux through this module is Jμ=yμ(S1-S2/K12), where

(1) yμ=F(xμ)

is the effective reaction rate constant of module μ (Figure 1B). Importantly, yμ depends only on the rate matrix xμ, but not on any other rate constants (see Corollary 2 in Materials and methods), and it can be recursively computed for any module, as described in Materials and methods. In other words, metabolic network 𝒩 can be coarse-grained by replacing module μ at steady state with a single first-order biochemical reaction with rate yμ. Importantly, such coarse-graining does not alter the dynamics of any metabolites outside of module μ (see Proposition 1 in Materials and methods). This statement is the biochemical analog of the star-mesh transformation (and its generalization, Kron reduction, Rao et al., 2014) well known in the theory of electric circuits (Versfeld, 1970). The biological interpretation of these properties is that a module is somewhat isolated from the rest of the metabolic network. And vice versa, the larger network (i.e. the cell) ‘cares’ only about the total rate at which the I/O metabolites are interconverted by the module but ‘does not care’ about the details of how this conversion is enzymatically implemented. In this sense, the effective rate yμ quantifies the function of module μ (a macroscopic parameter) while the rates xμ describe the specific biochemical implementation of the module (microscopic parameters).

The effective rate constant yμ of module μ depends on the entire rate matrix xμ. In general, a single mutation may perturb several rate constants within a module, so that the entire shape of the function F may be important. Here, I focus on a special case when each mutation perturbs one reaction (real or effective) within a module, while all others remain constant. To examine epistasis between mutations, I will also consider two different mutations that perturb two separate reactions within a module. In these special cases, we do not need to know the entire function F. We only need to know how module’s effective rate constant yμ depends on the one or two rate constants of the perturbed reactions. When yμ is considered as a function of the rate constant ξ of one reaction, I write

(2) yμ=f1(ξ),

and when yμ is considered as a function of the rate constants ξ and η of two reactions, I write

(3) yμ=f2(ξ,η).

The rate constants of all other reactions within module μ play a role of parameters in functions f1 and f2.

Consider now a network 𝒩 that has a hierarchical structure, such that there is a series of nested modules μν, in the sense that AμAν (Figure 1A). Since any module at steady state can be replaced with an effective first-order biochemical reaction, there exists a hierarchy of quantitative metabolic phenotypes yμ,yν, (Figure 1B,C). These phenotypes are of course functionally related to each other. Specifically, because ν is a ‘higher-level’ module (in the sense that it contains a ‘lower-level’ module μ), the matrix xν can be decomposed into two submatrices xμ and xνμ where the latter is the matrix of rate constants of reactions that belong to module ν but not to module μ. Since replacing the lower-level module μ with an effective reaction with rate constant yμ does not alter the dynamics of metabolites outside of μ, yν must depend on all elements of xμ only through yμ, that is,

(4) yν=f1(yμ),

where rates xνμ act as parameters of function f1. Thus, in the hierarchy of metabolic phenotypes yμ,yν,, a phenotype at each subsequent level depends on the phenotype at the preceding level according to Equation 4, and the lowest level phenotype yμ depends on the actual rate constants accroding to Equation 1. This hierarchy of functionally nested phenotypes is conceptually similar to the hierarchical ‘ontotype’ representation of genomic data proposed recently by Yu et al., 2016.

Quantification of epistasis

Consider a mutation A that perturbs only one rate constant xij, such that the wildtype value xij0 changes to xijA. This mutation can be quantified at the microscopic level by its relative effect δAxij=xijA/xij0-1. If the reaction between metabolites i and j belongs to nested modules μ,ν,, then mutation A may impact the functions of these modules, which can be quantified by the relative effects δAyμ, δAyν, etc. at each level of the hierarchy.

Consider now another mutation B that only perturbs the rate constant xk of another reaction. Since mutations A and B perturb distinct enzymes, they by definition do not genetically interact at the microscopic level. However, if both perturbed reactions belong to the metabolic module μ (and, as a consequence, to all higher-level modules which contain μ), they may interact at the level of the function of this module, in the sense that the effect of mutation B on the effective rate yμ may depend on whether mutation A is present or not. Such epistasis between mutations A and B can be quantified at the level μ of the metabolic hierarchy by a number of various epistasis coefficients (Wagner et al., 1998; Hansen and Wagner, 2001; Mani et al., 2008). I will quantify it with the epistasis coefficient

(5) εAByμ=δAByμ-δAyμ-δByμ2δAyμδByμ,

where δAByμ denotes the effect of the combination of mutations A and B on phenotype yμ relative to the wildtype. Since I only consider two mutations A and B, I will write εyμ instead of εAByμ to simplify notations. Note that other epistasis coefficients can always be computed from εyμ, δAyμ and δByμ, if necessary. Expressions for epistasis coefficients at other levels of the metabolic hierarchy are analogous.

Results

The central goal of this paper is to understand the patterns of epistasis between mutations that affect reaction rates in the hieararchical metabolic network described above. Specifically, I am interested in two questions. (1) Given that two mutations A and B have an epistasis coefficient εyμ at a lower level μ of the metabolic hierarchy, what can we say about their epistasis coefficient εyν at a higher level ν of the hierarchy? In other words, how does epistasis propagate through the metabolic hierarchy? (2) If mutation A only perturbs the activity xij of one enzyme and mutation B only perturbs the activity xk of another enzyme that belongs to the same module μ, then what values of εyμ can we expect to observe based on the topological relationship between the two perturbed reactions within module μ? In other words, what kinds of epistasis emerge in a metabolic network?

Propagation of epistasis through the hierarchy of metabolic phenotypes

Assuming that the effects of both individual mutations and their combined effect at the lower-level μ are small, it follows from Equation 4 and Equation 5 that

(6) εyν=εyμC+H2C2+o(1),

where C=f1yμ/yν and H=f1′′yμ2/yν are the first- and second-order control coefficients of the lower-level module μ with respect to the flux through the higher-level module ν and o(1) denotes all terms that vanish as the effects of mutations tend to zero (see Materials and methods for details). Note that Equation 6 is a special case of a more general Equation 49 which describes the case when mutations affect multiple enzymes. Equation 6 defines a linear map ϕ with slope 1/C and a fixed point ε¯=-H(2C(1-C))-1, which both depend on the topology of the higher-level module ν and the rate constants xνμ.

To gain some intuition for how the map ϕ governs the propagation of epistasis from a lower level μ to a higher level ν, suppose that module ν is a linear metabolic pathway. In this case, it is intuitively clear that function f1 is monotonically increasing (i.e. the higher yμ, the more flux can pass through the linear pathway ν) and concave (i.e. as yμ grows, other reactions in ν become increasingly more limiting, such that further gains in yμ yield smaller gains in yν). Indeed, it is easy to show that C=(1+αyμ)-1>0 and H=-2αyμ(1+αyμ)-2<0, where α is a positive constant that depends on other reactions in the pathway (see Materials and methods for details). It then immediately follows that any zero or negative epistasis εyμ that already arose at the lower level would propagate to negative epistasis εyν at the level of the linear pathway ν. Moreover, since C<1, the fixed point of the map in Equation 6 is unstable. Therefore, if epistasis εyμ was already sufficiently large at the lower level, it would induce even larger positive epistasis εyν at the level of the linear pathway ν. In fact, when module ν is a linear pathway, ε¯=1, so that εyν>1 whenever εyμ>1.

The first result of this paper is the following theorem, which shows that the same rules of propagation of epistasis hold not only for a linear pathway but for any module (Figure 2).

Propagation of epistasis.

Properties of Equation 6 that maps lower-level epistasis εyμ onto higher-level epistasis εyν. Slope 1/C and fixed point ε¯ depend on the topology and the rate constants of the higher-level module ν, but they are bounded, as shown. Thus, the fixed point ε¯ of this map lies between 0 and 1 and is always unstable (open circle).

Theorem 1

For any module ν,

(7) 0C1

and

(8) 0ε¯1.

The proof of Theorem 1 is given in Materials and methods. Its main idea is the following. The functional form of f1 in Equation 4 depends on the topology of module ν. Since the number of topologies of ν is infinite, we might a priori expect that there is also an infinite number of functional forms of f1. However, this is not the case. In fact, all higher-level modules that contain a lower-level module fall into three topological classes defined by the location of the lower-level module with respect to the I/O metabolites of the higher-level module (see Proposition 2 and Figure 7 in Materials and methods). To each topological class corresponds a parametric family of the function f1, so that there are only three such families. For each family, the values of C and H can be explicitly calculated, yielding the bounds in Equation 7 and Equation 8.

Equation 6 together with Equation 7 and Equation 8 show that the linear map ϕ from epistasis at a lower-level to epistasis at the higher-level has an unstable fixed point between 0 and 1 (Figure 2). This implies that negative epistasis at a lower level of the metabolic hierarchy necessarily induces negative epistasis of larger magnitude at the next level of the hierarchy, that is, εyνεyμ<0. Therefore, once negative epistasis emerges somewhere along the hierarchy, it will induce negative epistasis at all higher levels of the hierarchy, irrespectively of the topology or the kinetic parameters of the network.

Similarly, if epistasis at the lower level of the metabolic hierarchy is positive and strong, εyμ>1, it will induce even stronger positive epistasis at the next level of the hierarchy, that is, εyνεyμ>1. Therefore, once strong positive epistasis emerges somewhere in the metabolic hierarchy, it will induce strong positive epistasis of larger magnitude at all higher levels of the hierarchy, irrespectively of the topology or the kinetic parameters of the network. If positive epistasis at a lower level of the hierarchy is weak, 0<εyμ<1, it could induce either negative, weak positive or strong positive epistasis at the higher level of the hierarchy, depending on the precise value of εyμ, the topology of the higher-level module ν and the microscopic rate constants xνμ.

In summary, there are three regimes of how epistasis propagates through a hierarchical metabolic network. Negative and strong positive epistasis propagate robustly irrespectively of the topology and kinetic parameters of the metabolic network, whereas the propagation of weakly positive epistasis depends on these details. The strongest qualitative prediction that follows from Theorem 1 is that negative epistasis for a lower-level phenotype cannot turn into positive epistasis for a higher-level phenotype, but the converse is possible.

Emergence of epistasis between mutations affecting different enzymes

Which of the three regimes described above can emerge in metabolic networks and under what circumstances? In other words, if two mutations affect the same module, are there any constraints on epistasis that might arise at the level of the effective rate constant of this module? To address this question, I consider two mutations A and B that affect the rate constants of different single reactions within a given module.

Consider a relatively simple module ν shown in Figure 1A and two mutations A and B that affect the reactions, as shown in Figure 3A. I will now show that the epistasis coefficient εyν can take values in all three domains described above, depending on the biochemical details of this module. Using the recursive procedure for evaluating yμ described in Materials and methods, it is straightforward to obtain an analytical expression for yν as a function of the rate matrix xν, from which εyμ can also be obtained (see Materials and methods for details). To demonstrate that εyμ can take values below 0, between 0 and 1, and above 1, it is convenient to keep all of the rate constants fixed except for the rate constant zx34 of a reaction that is not affected by mutations A or B, as shown in Figure 3A. Figure 3B then shows how the epistasis coefficient εyμ varies as a function of z for one particular choice of all other rate constants. When z is small, εyμ<0. As z increases, it becomes weakly positive (0<εyμ<1) and eventually strongly positive (εyμ>1). Thus, in my model, there are no fundamental constraints on the types of epistasis that can emerge between mutations.

Emergence of epistasis and its dependence on the topological relationship between the reactions affected by mutations.

(A) An example of a simple module ν (same as in Figure 1A) where negative, weak positive and strong positive epistasis can emerge between two mutations A and B. (B) Epistasis between mutations A and B at the level of module ν depicted in (A) as a function of the rate constant z of a third reaction. The values of other parameters of the network are given in Materials and Methods. (C) An example of a simple module where reactions affected by mutations are strictly parallel. In such cases, epistasis for the effective rate constant yν is non-positive. Dashed blue lines highlight paths that connect the I/O metabolites and each contain only one of the affected reactions. (D) An example of a simple module where reactions affected by mutations are strictly serial. In such cases, epistasis for the effective rate constant yν is equal to or greater than 1 (i.e. strongly positive). Dashed blue line highlights a path that connects the I/O metabolites and contains both affected reactions.

This simple example also reveals that not only the value but also the sign of epistasis generically depend on the rates of other reactions in the network, such that other mutations or physiological changes in enzyme expression levels can modulate epistasis sign and strength. In other words, ‘higher-order’ and ‘environmental’ epistasis are generic features of metabolic networks.

Upon closer examination, the toy example in Figure 3 also suggests that the sign of εyν may depend predictably on the topological relationship between the affected reactions. When z=0, the two reactions affected by mutations are parallel, and epistasis is negative. When z is very large, most of the flux between the I/O metabolites passes through z such that the two reactions affected by mutations become effectively serial, and epistasis is strongly positive. Other toy models show consistent results: epistasis between mutations affecting different reactions in a linear pathway is always positive and epistasis between mutations affecting parallel reactions is negative (see Materials and methods for details). These observation suggest an interesting conjecture. Do mutations affecting parallel reactions always exhibit negative epistasis and do mutations affecting serial reactions always exhibit positive epistasis? In fact, such relationship between sign of epistasis and topology has been previously suggested in the literature (e.g. Dixon et al., 2009; Lehner, 2011).

To formalize and mathematically prove this hypothesis, I first define two reactions as parallel within a given module if there exist at least two distinct simple (i.e. non-self-intersecting) paths that connect the I/O metabolites, such that each path lies within the module and contains only one of the two focal reactions. Analogously, two reactions are serial within a given module if there exists at least one simple path that connects the I/O metabolites, lies within the module and contains both focal reactions.

According to these definitions, two reactions can be simultaneously parallel and serial, as, for example, the reactions affected by mutations A and B in Figure 3A. I call such reaction pairs serial-parallel. I define two reactions to be strictly parallel if they are parallel but not serial (Figure 3C) and I define two reactions to be strictly serial if they are serial but not parallel (Figure 3D). Thus, each pair of reactions within a module can be classified as either strictly parallel, strictly serial or serial-parallel.

The second result of this paper is the following theorem.

Theorem 2

Let ξ and η be the rate constants of two different reactions in module μ. Suppose that mutation A perturbs only one of these reactions by δAξ and mutation B perturbs only the other reaction by δBη. In the limit δAξ0 and δBη0, the following statements are true. If the affected reactions are strictly parallel then εyμ0. If the affected reactions are strictly serial, then εyμ1.

The detailed proof of this theorem is given in Materials and methods. Its key ideas and the logic are the following. It follows from Equation 3 and Equation 5 that

(9) εyμ=Hξη2CξCη+o(1),

where Cξ=f2ξξyμ, Cη=f2ηηyμ, Hξη=2f2ξηξηyμ are the first- and second-order control coefficients of the affected reactions with respect to the flux through module μ and o(1) denotes terms that vanish when δAξ and δBη approach zero (see Materials and methods for details). Note that Equation 9 was previously derived by Chiu et al., 2012.

To compute the epistasis coefficient εyμ for an arbitrary module μ, we need to know the first and second derivatives of function f2. Analogous to function f1, there is a finite number of parametric families to which f2 can belong. Specifically, all modules fall into nine topological classes with respect to the locations of the affected reactions within the module (see Figure 8), and each of these topologies defines a parametric family of function f2 (see Proposition 3 and its Corollary 3 in Materials and methods). Most of these topological classes are broad and contain modules where the affected reactions are strictly parallel, those where they are strictly serial as well as those where they are serial-parallel. And it is easy to show that not all members of each topological class have the same sign of εyμ. However, modules from the same topological class where the affected reactions are strictly parallel or strictly serial fall into a finite number of topological sub-classes (see Figure 10 through Figure 14, Table 2 and Table 3). Overall, there are only 17 distinct topologies where the affected reactions are strictly parallel (Table 2), which define 17 parametric sub-families of function f2. For all members of these sub-families, Equation 9 yields εyμ0 (see Proposition 7 in Materials and methods). Similarly, there are only 11 distinct topologies where the affected reactions are strictly serial (Table 3), which define 11 parametric sub-families of function f2. For all members of these sub-families, Equation 9 yields εyμ1 (see Proposition 8 in Materials and methods).

The results of Theorem 1 and Theorem 2 together imply that the topological relationship at the microscopic level between two reactions affected by mutations constrains the values of their epistasis coefficient at all higher phenotypic levels. Specifically, if negative epistasis is detected at any phenotypic level, the affected reactions cannot be strictly serial. And conversely, if strong positive epistasis is detected at any phenotypic level, the affected reactions cannot be strictly parallel. In this model, weak positive epistasis in the absence of any additional information does not imply any specific topological relationship between the affected reactions.

Sensitivity of results with respect to the magnitude of mutational effects

Both Theorem 1 and Theorem 2 strictly hold only when the effects of both mutations are infinitesimal. Next, I investigate how these results might change when the mutational effects are finite.

Propagation of epistasis between mutations with finite effect sizes

As mentioned above and discussed in detail in Materials and methods, all higher-level modules that contain a lower-level module fall into three topological classes, which I label b, io and i, depending on the location of the lower-level module within the higher- level module (see Figure 7). The topological class specifies the parametric family of the function f1 which maps the effective rate constant yμ onto the effective rate constant yν (see Equation 4). For all modules from the topological class b, function f1 is linear (see Equation 29), which implies that the results of Theorem 1 hold exactly even when the effects of mutations are finite. For modules from the topological classes io and i, function f1 is hyperbolic (see Equation 30 and Equation 31), so that the results of Theorem 1 may not hold when the effects of mutations are finite. To test the validity of Theorem 1 in these cases, I calculated the non-linear function ϕ that maps the epistasis coefficient εyμ onto the epistasis coefficient εyν for 1000 randomly generated modules from each of the two topological classes and for mutations that increase or decrease the effective rate constant of the lower-level module yμ by up to 50% (see Materials and methods for details).

The validity of Theorem 1 depended on the sign of mutational effects. When at least one of the two mutations had a negative effect on yμ, map ϕ had the same properties as described in Theorem 1, even for mutations with large effect, that is, it had a fixed point ε¯ in the interval [0,1] and this fixed point was unstable. When the effects of both mutations on yμ were positive and small, these results also held in about 82% of sampled modules (see Figure 4A, Figure 4—figure supplement 1, Figure 4—figure supplement 2). In the remaining ∼18% of sampled modules, the fixed point ε¯ shifted slightly above 1. As the magnitude of mutational effects increased, the fraction of sampled modules with ε¯>1 grew, reaching 42% when both mutations increased yμ by 50%. In most of these cases, ε¯ remained below 2, and I found only one module with ε¯>4 (Figure 4A, Figure 4—figure supplement 1, Figure 4—figure supplement 2). Whenever the fixed point existed, it was unstable, with the exception of 12 modules for which ϕ was very close to the identity map. For 289 modules (14.5%), the fixed point disappeared when both mutations increased yμ by 50%. In all these cases, εyν<εyμ, indicating that even large positive epistasis may decline as it propagates through the metabolic hierarchy when the effects of mutations are finite.

Figure 4 with 3 supplements see all
Sensitivity of results of Theorem 1 and Theorem 2 with respect to the magnitude of mutational effects.

(A) Distribution of the position of the fixed point ε¯ of the function ϕ that maps lower-level epistasis εyμ onto higher-level epistasis εyν in modules with random parameters and for mutations with positive effects on yμ (see text and Materials and methods for details). All cases are shown in Figure 4—figure supplement 1 and Figure 4—figure supplement 2. The effect size of both mutations is indicated on each panel. ‘No f.p'. indicates that no fixed point exists. (B) Fraction of sampled modules (averaged across generating topologies) where mutations affect strictly serial reactions but the epistasis coefficient is less than 1, contrary to the statement of Theorem 2 (see text and Materials and methods for details). All cases stratified by generating topology are shown in Figure 4—figure supplement 3.

Emergence of epistasis between mutations with finite effect sizes

As mentioned above and discussed in detail in Materials and methods, modules where the reactions affected by mutations are strictly parallel fall into 17 topological classes (see Table 2) and modules where the reactions affected by mutations are strictly serial fall into 11 topological classes (see Table 3). The topological class specifies the parametric family of the function f2 which maps the rate constants ξ and η of the affected reactions onto the effective rate constant yμ. To test how well Theorem 2 holds when the effects of mutations are finite, I calculated εyμ for randomly generated modules from these topological classes and for mutations increasing or decreasing ξ and η by up to 50% (see Materials and methods for details).

The validity of Theorem 2 depended most strongly on the topological relationship between the reaction affected by mutations. Whenever the affected reactions were strictly parallel, the epistasis coefficient at the level of module μ was always less than or equal to zero, even when mutations perturbed the rate constants by as much as 50%, consistent with Theorem 2. This was also true for strictly serial reactions, as long as both mutations had positive effects. When the affected reaction were strictly serial and at least one of the mutations had a negative effect, the epistasis coefficient was always positive, but in some cases it was less than 1 (see Figure 4B, Figure 4—figure supplement 3), in disagreement with Theorem 2. This indicates that when the effects of mutations are not infinitesimal, even mutations that affect strictly serial reactions can potentially produce negative epistasis for higher-level phenotypes.

Taken together, these results suggest that both Theorem 1 and Theorem 2 extend reasonably well, but not perfectly, to mutations with finite effect sizes. The domains of validity of both theorems appear to depend on the sign of mutational effects. The way in which the theorems break down as their assumptions are violated appears to be stereotypical: when the mutational effects increase, more types of mutations produce weak epistasis, and the bias toward negative epistasis increases during propagation from lower to higher levels of the metabolic hierarchy.

Beyond first-order kinetics: epistasis in a kinetic model of glycolysis

The results of previous sections revealed a relationship between network topology and the ensuing epistasis coefficients in an analytically tractable model. However, the assumptions of this model are most certainly violated in many realistic situations. It is therefore important to know whether the same or similar rules of epistasis emergence and propagation hold beyond the scope of this model. I address this question here by analyzing a computational kinetic model of glycolysis developed by Chassagnole et al., 2002. This model keeps track of the concentrations of 17 metabolites, reactions between which are catalyzed by 18 enzymes (Figure 5A and Figure 5—figure supplement 1; see Materials and methods for details). This model falls far outside of the analytical framework introduced in this paper: some reactions are second-order, reaction kinetics are non-linear, and in several cases the reaction rates are modulated by other metabolites (Chassagnole et al., 2002).

Figure 5 with 3 supplements see all
Epistasis in a kinetic model of Escherichia coli glycolysis.

(A) Simplified schematic of the model (see Figure 5—figure supplement 1 for details). Different shades of gray in the background highlight four modules as indicated (see text). Light blue circles represent metabolites. Reactions are shown as lines with dark gray boxes. The enzymes catalyzing reactions whose control coefficients with respect to the flux through the module are positive are named; other enzyme names are ommitted for clarity (see Table 5 and Table 6 for abbreviations). Three reactions, catalyzed by PGI, PFK, PGDH, for which the epistasis coefficients are shown in panel B are highlighted in dark blue, red, and orange, respectively. (B) Epistasis coefficients for flux through each module between mutations perturbing the respective reactions, computed at steady state (see text and Materials and methods for details). Reactions catalyzed by PGI and PGDH are strictly parallel (path g6p-f6p-fdp-gap contains only PGI, path g6p-6pg-ribu5p-gap contains only PGDH and there is no simple path in UGPP between g6p and gap that contains both PGI and PGDH). Reactions catalyzed by PGI and PFK are serial-parallel (path g6p-f6p-fdp-gap contains both reactions, path g6p-f6p-gap contains only PGI, path g6p-6pg-ribu5p-f6p-fdp-gap contains only PFK). Reactions catalyzed by PFK and PGDH are also serial-parallel (path g6p-6pg-ribu5p-f6p-fdp-gap contains both reactions, path g6p-f6p-fdp-gap contains only PFK, path g6p-6pg-ribu5p-gap contains only PGDH).

Testing the predictions of the analytical theory in this computational model faces two complications. First, in a non-linear model, modules are no longer fully characterized by their effective rate constants, even at steady state. Instead, each module is described by the flux between its I/O metabolites which non-linearly depends on the concentrations of these metabolites. Consequently, the effects of mutations and epistasis coefficients also become functions of the I/O metabolite concentrations. An epistasis coefficient at the level of module ν can still be evaluated according to Equation 5, with yν now representing the flux through module ν evaluated at a particular concentration of the I/O metabolites. For simplicity, I computationally find the steady state of the full glycolysis network and evaluate the epistasis coefficients only at this steady state, that is, for each module, I keep the concentrations of the I/O metabolites fixed at their steady-state values for the full network (see Materials and methods for details).

The second complication is that some control coefficients are so small that they fall below the threshold of numerical precision. Perturbing such reactions has no detectable effect on flux (Figure 5—figure supplement 2). In the analysis that follows, I ignore such reactions because the epistasis coefficient defined by Equation 5 can only be computed for mutations with non-zero effects on flux. In addition, the control coefficients of some reactions are negative, which implies that an increase in the rate of such reaction decreases the flux through the module (Figure 5—figure supplement 2). I also ignore such reactions because there is no analog for them in the analytical theory presented above. After excluding seven reactions for these reasons, I examine epistasis in 55 pairs of mutations that affect the remaining 11 reactions.

The glycolysis network shown in Figure 5A (see also Figure 5—figure supplement 1) can be naturally partitioned into four modules which I name ‘LG’ (lower glycolysis), ‘UGPP’ (upper glycolysis and pentose phosphate), ‘GPP’ (glycolysis and pentose phosophate), and ‘FULL’. Modules LG and UGPP are non-overlappng and both of them are nested in module GPP which in turn is nested in the FULL module. Thus, at least for some reaction pairs it is possible to calculate epistasis coefficients at three levels of metabolic hierarchy. There are three such pairs, and the results for them are shown in Figure 5B. Epistasis for the remaining pairs of reactions can be evaluated only at one or two levels of the hierarchy because these reactions belong to different modules at the lowest levels or because their individual effects are too small. The results for all reaction pairs are shown in Figure 5—figure supplement 3.

The strongest qualitative prediction of the analytical theory described above is that negative epistasis for a lower-level phenotype cannot turn into positive epistasis for a higher-level phenotype, while the converse is possible. Figure 5B and Figure 5—figure supplement 3 show that the data are consistent with this prediction. Another prediction is that epistasis between strictly parallel reactions should be negative. There is only one pair of reactions that are strictly parallel, those catalyzed by glucose-6-phosphate isomerase (PGI) and 6-phosphogluconate dehydrogenase (PGDH), and indeed the epistasis coefficients between mutations affecting these reactions are negative at all levels of the hierarchy (Figure 5B). Finally, the analytical theory predicts that mutations affecting strictly serial reactions should exhibit strong positive epistasis. There are 36 reaction pairs that are strictly serial. Epistasis is positive between mutations in 33 of them, and it is strongly positive in 17 of them (Figure 5—figure supplement 3). Three pairs of strictly serial reactions (those where one reaction is catalyzed by PK and the other is catalyzed by PGI, PGDH, or PFK) exhibit negative epistasis (Figure 5—figure supplement 3). These results suggest that, although one may not be able to naively extrapolate the rules of emergence and propagation of epistasis derived in the simple analytical model to more complex networks, some generalized versions of these rules may nevertheless hold more broadly.

Discussion

Genetic interactions are a powerful tool in genetics, and they play an important role in evolution. Yet, how epistasis emerges from the molecular architecture of the cell and how it propagates to higher-level phenotypes, such as fitness, remains largely unknown. Several recent studies made a statistical argument that the structure of the fitness landscape (and, as a consequence, the epistatic interactions between mutations at the level of fitness) may be largely independent of the underlying molecular architecture of the organism (Martin, 2014; Lyons et al., 2020; Reddy and Desai, 2020). If mutations are typically highly pleiotropic (i.e. affect many independent phenotypes relevant for fitness) or are engaged in a large number of idiosyncratic epistatic interactions with other mutations in the genome, the resulting fitness landscapes converge to certain limiting shapes, such as the Fisher’s geometric model (Martin, 2014; Tenaillon, 2014). To what extent these arguments indeed apply in practice is unclear. But if they do, most genetic interactions detected at the fitness level may be uninformative about the architecture of the underlying biological networks.

In this paper, I took a ‘mechanistic’ approach, which is in a sense orthogonal to the statistical one. In my model of a hierarchical metabolic network, mutations are highly pleiotropic (a mutation in any enzyme affects all the fluxes in the module) and highly epistatic (a mutation in any enzyme interacts with mutations in any other enzyme). Yet, these pleiotropic and epistatic effects appear to be sufficiently structured that some information about the topology of the network is preserved through all levels of the hierarchy. Indeed, the emergence and propagation of epistasis follow two simple rules in my model. First, once epistasis emerges at some level of the hierarchy, its propagation through the higher levels of the hierarchy depends weakly on the details of the network. Specifically, negative epistasis at a lower level induces negative epistasis at all higher levels and strong positive epistasis induces strong positive epistasis at all higher levels, irrespectively of the topology or the kinetic parameters of the network. Second, what type of epistasis emerges in the first place depends on the topological relationship between the reactions affected by mutations. In particular, negative epistasis emerges between mutations that affect strictly parallel reactions and positive epistasis emerges between mutations that affect strictly serial reactions. Insofar as my model is relevant to nature, the key conclusion from it is that epistasis at high-level phenotypes carries some, albeit incomplete, information about the underlying topological relationship between the affected reactions.

These results have implications for the interpretation of empirically measured epistasis coefficients. It is often assumed that a positive epistasis coefficient between mutations that affect distinct genes signals that their gene products act in some sense serially, whereas a negative epistasis coefficient is a signal of genetic redundency, that is, a parallel relationship between gene products (Dixon et al., 2009). My results suggest that this reasoning is generally correct, but that the relationship between epistasis and topology is more nuanced. In particular, the sign of the epistasis coefficient in my model constrains but does not uniquely specify the topological relationship, such that a negative epistasis coefficient implies that the affected reactions are not strictly serial (but may or may not be strictly parallel) and an epistasis coefficient exceeding unity excludes a strictly parallel relationship (but does not necessarily imply a strictly serial relationship). My model suggests that one should also be careful with inferences going in the other direction, that is, extrapolating the patterns of epistasis measured at the biochemical level to those at the level of fitness. For example, if one wishes to infer the past evolutionary trajectory of an enzyme and finds two amino acid changes that exhibit a positive interaction at the level of enzymatic activity, it does not automatically imply that these mutations will exhibit a positive interaction at the level of fitness.

The strongest results presented here rely on several assumptions. I proved Theorem 1 and Theorem 2 in the limit of vanishingly small mutational effects. Some results of the metabolic control analysis, notably the summation theorem, are sensitive to this assumption (Bagheri-Chaichian et al., 2003; Bagheri and Wagner, 2004). To test the sensitivity of my analytical results with respect to this assumption, I used numerical simulations of networks with randomly sampled kinetic parameters and found that the results hold reasonably well when the effects of mutations are not infinitesimal.

The most restrictive assumption in the present work is that of first-order kinetics. Networks with only first-order kinetics clearly fail to capture some biologically important phenomena, such as sign epistasis (Weinreich et al., 2005; Chou et al., 2014; Ewald et al., 2017; Kemble et al., 2020). I discuss possible ways to relax this assumption below. But at present, a major question remains whether the rules of epistasis and propagation described here hold for realistic biological networks and whether they can be directly used to interpret empirical epistasis coefficients. My analysis of a fairly realistic computational model of glycolysis cautions against overinterpreting empirical epistasis coefficients using the rules derived here. But it also suggests that more general rules of propagation and emergence of epistasis may be found for more realistic networks. Thus, the simple rules derived here should probably be thought of as null expectations.

Relaxing the first-order kinetics assumption is analytically challenging because it is critical for replacing a module with a single effective reaction without altering the dynamics of the rest of the network. Although such lossless replacement is almost certainly not possible in networks with more complex kinetics, advanced network coarse-graining techniques may offer a promising way forward (Rao et al., 2014). Flux balance analysis (FBA) is an alternative approach (Orth et al., 2010). FBA is appealing because it entirely removes the dependence of the model on reaction kinetics. However, this comes at a substantial cost. In FBA models, fitness and other high-level phenotypes become independent of the internal kinetic parameters, which is clearly unrealistic. Nevertheless, FBA is often very good at capturing the effects of mutations that change the topology of metabolic networks, such as reaction additions and deletions (reviewed in Gu et al., 2019). At the same time, there is no natural way within FBA to model mutations that perturb reaction kinetics (He et al., 2010; Alzoubi et al., 2019). In short, FBA and my approach are complementary (see Appendix 5 for a more detailed discussion).

Generic properties of epistasis in biological systems

Simple models help us identify generic phenomena—those that are shared by a large class of systems—which should inform our ‘null’ expectations in empirical studies. Deviations from such null in a given system under examination inform us about potentially interesting peculiarities of this system. The model presented here suggests several generic features of epistasis between genome-wide mutations.

Epistasis has two contributions

My analysis shows that the value of an epistasis coefficient measured for a higher level phenotype is a result of two contributions (Domingo et al., 2019), propagation and emergence, which correspond to two terms in Equation 6 (or the more general Equation 49). The first term, propagation, shows that if two mutations exhibit epistasis for a lower-level phenotype they also generally exhibit epistasis for a higher-level phenotype. The second contribution comes from the fact that lower-level phenotypes map onto higher-level phenotypes via non-linear functions. This is true even in a simple model with linear kinetics considered here. As a result, even if two mutations exhibit no epistasis at the lower-level phenotype, epistasis must emerge for the higher-level phenotype, as previously pointed out by multiple authors (e.g. Kacser and Burns, 1981; DePristo et al., 2005; Martin et al., 2007; Chiu et al., 2012; Otwinowski et al., 2018; Domingo et al., 2019; Husain and Murugan, 2020).

Epistasis depends on the genetic background and environment

My analysis shows that the value of an epistasis coefficient for a particular pair of mutations is in large part determined by the topological relationship between reactions affected by them. Since the topology of the metabolic network itself depends on the genotype (which genes are present in the genome) and on the environment (which enzymes are active or not), the topological relationship between two specific reactions might change if, for example, a third mutation knocks out another enzyme or if an enzyme is up- or down-regulated due to an environmental change (see Figure 3). Thus, we should generically expect epistasis between mutations to depend on the environment and on the presence or absence of other mutations in the genome. In other words, G×G×G interactions (higher-oder epistasis) and G×G×E interactions (environmental epistasis) should be common (Snitkin and Segrè, 2011; Flynn et al., 2013; Lindsey et al., 2013; Taylor and Ehrenreich, 2015; Sailer and Harms, 2017a). This fact complicates the interpretation of inter-gene epistasis since mutations in the same pair of genes can exhibit qualitatively different genetic interactions in different strains, organisms and environments, as has been observed (St Onge et al., 2007; Musso et al., 2008; Tischler et al., 2008; Dowell et al., 2010; Heigwer et al., 2018; Li et al., 2019). However, the situation may not be hopeless because the topological relationship between two reactions cannot change arbitrarily after addition or removal of a single reaction. For example, if two reactions are strictly parallel, removing a third reaction does not alter their relationship (see Proposition 5). Thus, comparing matrices of epistasis coefficients measured in different environments or genetic backgrounds could inform us about how the organism rewires its metabolic network in response to these perturbations (St Onge et al., 2007; Musso et al., 2008; Heigwer et al., 2018; Li et al., 2019).

Skew in the distribution of epistasis coefficients

Studies that measure epistasis for fitness-related phenotypes among genome-wide mutations usually find both positive and negative epistases, but the preponderance of positive and negative epistasis varies. Some authors reported a skew toward positive interactions among deleterious mutations (Jasnos and Korona, 2007; He et al., 2010; Johnson et al., 2019), whereas others reported a skew toward negative interactions (Szappanos et al., 2011; Costanzo et al., 2016). Beneficial mutations appear to predominantly exhibit negative epistasis, also known as ‘diminishing returns’ epistasis (e.g. Martin et al., 2007; Khan et al., 2011; Chou et al., 2011; Kryazhimskiy et al., 2014; Schoustra et al., 2016). The reasons for these patterns are currently unclear. Several recent theoretical papers offer possible statistical explanations for them (Martin, 2014; Lyons et al., 2020; Reddy and Desai, 2020). On the other hand, mechanistic predictions for the distribution of epistasis coefficients are not yet available (but see Sanjuán and Nebot, 2008; Macía et al., 2012; Chiu et al., 2012). The present work does not directly address this problem either, but it provides some additional clues.

First, my model shows that the sign of epistasis at least to some extent reflects the topology of the network. Thus, the distribution of epistasis coefficients at high-level phenotypes in real organisms should ultimately depend on the preponderance of different topological relationships between the edges in biological networks. It then seems a priori unlikely that positive and negative interactions would be exactly balanced. Thus, we should expect the distribution of epistasis coefficients to be skewed in one or another direction.

The second observation is that in the metabolic model considered here a positive epistasis coefficient at one level of the hierarchy can turn into a negative one at a higher level, but the reverse is not possible. This bias toward negative epistasis at higher-level phenotypes appears to be even stronger in networks with saturating kinetics (Figure 5 and Figure 5—figure supplement 3).

The third observation is that epistasis among beneficial and deleterious mutations affecting metabolic genes should be identical at the level where they arise, provided that beneficial and deleterious mutations are identically distributed among metabolic reactions. Thus, a stronger skew toward negative epistasis among beneficial mutations at the level of fitness could arise in my model for two mutually non-exclusive reasons. One possibility is that beneficial mutations tend to affect certain special subsets of genes, those that predominantly give rise to negative epistasis. For example, beneficial mutations may for some reason predominantly arise in enzymes that catalyze strictly parallel reactions. Another possibility is that when epistasis between beneficial mutations propagates through the metabolic hierarchy it tends to exhibit a stronger negative bias compared to epistasis between deleterious mutations. Indeed, this phenomenon arises in my model among mutations with large effect (see Figure 4A, Figure 4—figure supplement 1 and Figure 4—figure supplement 2).

Epistasis is generic

Perhaps the most important—and also the most intuitive—conclusion of this work is that we should expect epistasis for high-level phenotypes, such as fitness, to be extremely common. Consider first a unicellular organism growing exponentially. Its fitness is fully determined by its growth rate, which can be thought of as the rate constant of an effective biochemical reaction that converts external nutrients into cells (see Appendix 6 for a simple mathematical model of this statement). In other words, growth rate is the most coarse-grained description of a metabolic network and, as such, it depends on the rate constants of all underlying biochemical reactions. Many previous studies have shown that within-protein epistasis is extremely common (e.g. Lunzer et al., 2005; DePristo et al., 2005; Sailer and Harms, 2017b; Husain and Murugan, 2020). Present work shows that, once epistasis arises at the level of protein activity, it will propagate all the way up the metabolic hierarchy and will manifest itself as epistasis for growth rate. It also suggests that growth rate is a generically non-linear function of the rate constants of the underlying biochemical reactions, such that all mutations that affect growth rate individually would also exhibit pairwise epistasis for growth rate with each other (Kacser and Burns, 1981; DePristo et al., 2005; Martin et al., 2007; Chiu et al., 2012; Otwinowski et al., 2018; Husain and Murugan, 2020).

In more complex organisms and/or in certain variable environments, it may be possible to decompose fitness into multiplicative or additive components, for example, plant’s fitness may be equal to the product of the number of seeds it produces and their germination probability, as pointed out by Chou et al., 2011. Then, mutations that affect different components of fitness would exhibit no epistasis. However, such situations should be considered exceptional, as they require fitness to be decomposable and mutations to be non-pleiotropic.

If epistasis is in fact generic for high-level phenotypes, why do we not observe it more frequently? For example, a recent study that tested almost all pairs of gene knock-out mutations in yeast found genetic interactions for fitness for only about 4% of them (Costanzo et al., 2016). One possibility is that many pairs of mutations exhibit epistasis that is simply too small to detect with current methods. As the precision of fitness measurements improves, we would then expect the fraction of interacting gene pairs to grow. Another possibility is that systematic shifts in the distribution of estimated epistasis coefficients away from zero are taken by researchers as systematic errors rather than real phenomena, and are normalized out. Thus, some epistasis that would otherwise be detectable may be lost during data processing.

If epistasis is indeed as ubiquituous as the present analysis suggests, it would call into question how observations of inter-gene epistasis are interpreted. In particular, contrary to a common belief, a non-zero epistasis coefficient does not necessarily imply any specific functional relationship between the components affected by mutations beyond the fact that both components somehow contribute to the measured phenotype (Boyle et al., 2017). The focus of future research should then be not merely on documenting epistasis but on developing theory and methods for a robust inference of biological relationships from measured epistasis coefficients.

Materials and methods

Key ideas and logic of proofs of Theorems 1 and 2

Request a detailed protocol

Before proceeding to the detailed proofs of Theorem 1 and Theorem 2, I informally outline some key ideas and the basic logic.

The central object of the theory is a metabolic module. Modules have two key properties. First, a module is somewhat isolated from the rest of the metabolic network, in the sense that all metabolites inside it interact with only two metabolites outside, the I/O metabolites. The second property is that the metabolites within the module are sufficiently connected that each of them individually as well as any subset of them collectively can achieve a quasi-steady state (QSS), given the concentrations of the remaining metabolites. This property is proven in Proposition 1.

When some metabolites are at QSS, they can be effectively removed from the network and replaced with effective reactions among the remaining metabolites. In other words, one can ‘coarse-grain’ the network by removing metabolites. This approach is a standard biochemical-network reduction technique (Segel, 1988); for example, the Briggs-Haldane derivation of the Michaelis-Menten formula is based on this idea. In general, the resulting effective reactions have more complex (non-linear) kinetics than the original reactions. However, when the original reactions are first-order, the effective reactions are also first-order, that is, there is no increase in complexity. In Network coarse graining and an algorithm for evaluating the effective rate constant for an arbitrary module, I formally define the coarse-graining procedure (CGP) that eliminates one or multiple metabolites and replaces them with effective reactions.

CGP is an essential concept in my theory. I use it to compute the QSS concentrations for internal metabolites within a module (Corollary 1) and thereby prove Proposition 1, mentioned above. Since any module μ can achieve a QSS at any concentrations of its I/O metabolites and since any module has only two I/O metabolites, it can be replaced with a single effective reaction (Corollary 2). CGP provides a way to calculate the rate constant yμ of this reaction. In other words, the CGP is an algorithm for evaluating function F(xμ) in Equation 1 for any module μ.

CGP has an important property: its result does not depend on the order in which metabolites are eliminated. Therefore, in computing the effective rate constant of a module, we can choose any convinient way to eliminate its metabolites. Suppose that one module μ is nested within another module ν as in Figure 1A. A convenient way to compute the effective rate yν of the larger module is to first coarse-grain the smaller module μ, replacing it with an effective rate yμ, and then eliminate all the remaining metabolites in ν. Since effective rates after coarse-graining do not depend on the order of metabolite elimination, yν must depend on the rate constants xμ only indirectly, through yμ. In other words, all the information about the smaller module μ that is relevant for the performance of the larger module ν is contained in yμ. Therefore, if a mutation or mutations perturb only reactions inside of the smaller module μ, we only need to know their effects on the effective rate constant yμ to completely understand how they will perturb the performance of the larger module ν. Specifically, if we have two such mutations A and B, all the information about them is contained in three numbers, δAyμ, δByμ and εyμ. Theorem 1 then describes how epistasis at the level of module μ propagates to epistasis at the level of module ν.

The proof of Theorem 1 proceeds as follows. Let a be the effective reaction with rate constant yμ that represents module μ within the larger module ν, and consider yν as a function of yμ, as in Equation 4. To obtain f1(yμ), it is convenient to first eliminate all metabolites that do not participate in reaction a. No matter what the initial structure of module ν is, such coarse-graining will produce only one of three topologically distinct ‘minimal’ modules, which differ by the location of reaction a with respect to the I/O metabolites of module ν (Figure 7). This implies that the function f1 can belong to three parameteric families, where the parameters are the effective rate constants of reactions other than a in each of the minimal modules. This is the essence of Proposition 2. Then Theorem 1 can be easily proven by explicitly evaluating the first- and second-order control coefficients for each of the three functions and showing that the statements of the theorem hold for all of them, irrespectively of the function’s parameters.

Now consider two reactions a and b with rate constants ξ and η, and imagine the two mutations A and B that affect these reactions. To understand what value of εyμ will occur, we need to obtain yμ as a function of ξ and η (Equation 3). To do so, it is convenient to first eliminate all metabolites that do not participate in reactions a or b. No matter what the initial structure of module μ is, such coarse-graining will produce only one of nine topologically distinct minimal modules, which differ by the location of reactions a and b with respect to the I/O metabolites of module μ and each other (Figure 8). This implies that the function f2 can belong to nine parameteric families. This is the essence of Proposition 3 and Corollary 3.

How does the topological relationship between reactions a and b translate into epistasis? First, there are only three types of relationships between any pair of reactions in a module: strictly serial, strictly parallel, or serial-parallel (see Figure 3). Second, Proposition 4 and Corollary 4 show that coarse-graining does not alter the strict relationships, that is, if reactions a and b are strictly serial or strictly parallel before coarse-graining they will remain so after coarse-graining. This is important because it implies that to prove Theorem 2 we do not need to consider an infinitely large space of all modules but only a much smaller space of all minimal modules, that is, those that have only those metabolites that participate in the affected reactions a and b. Although the space of all minimal modules is still infinite, the space of their topologies is finite (see Figure 8). For some minimal topologies, the connection between the strictly serial or strictly parallel relationship and epistasis can be established very easily. For example, if reaction a and reaction b both share an I/O metabolite as a substrate (see topological class io,io,IO in Figure 8), then they are always strictly parallel, no matter what the rest of the module looks like. Evaluating the first- and second-order control coefficients for the function f2 that corresponds to this topological class reveals that εyμ0 for any parameter values of f2.

Unfortunately, most topological classes are too broad and include modules where reactions a and b are strictly serial as well as modules where they are strictly parallel or serial-parallel (e.g., class io,io,). Consequently, the sign of εyμ for such modules can change depending on the values of the rate constants. However, since the number of distinct minimal topologies is finite, it is possible to identify all minimal topologies where the reactions a and b are strictly serial or strictly parallel. These topological sub-classes define parametric sub-families of function f2, and we can explicitly calculate εyμ for all such functions. However, such brute-force approach is extremely cumbersome because the number of distinct minimal topologies is very large.

Fortunately, the following simple and intuitive fact greatly simplifies this problem. If two reactions are strictly serial or strictly parallel, this relationship does not change if a third reaction is removed from the module. This statement is the essence of Proposition 5. However, if the two reactions are serial-parallel, removal of a third reaction can change the relationship to a strictly serial or a strictly parallel one. As a consequence, there exist certain most connected ‘generating’ topologies where the relationship between the focal reactions is strictly parallel or strictly serial, and any other strictly serial minimal topology can be produced from at least one of the generating topologies by removal of reactions. This is the essence of Proposition 6. All generating topologies can be discovered by a simple algorithm provided in Appendix 3. They are listed in Table 2 and Table 3 and shown in Figure 10 through Figure 14. Each generating topology defines a parametric sub-family of function f2, and I explicitly evaluate the first- and second-order control coefficients for all these sub-families (see Proposition 7 and Proposition 8) which essentially completes the proof of Theorem 2.

Network coarse-graining

Notations and definitions

Request a detailed protocol

Here, I give a more precise definition of the model and introduce additional notations and definitions. As mentioned above, all reactions are first order and reversible. Thus, each reaction ij has one substrate iA and one product jA, and it is fully described by its rate constant xij. By definition, xii=0. I denote the set of all reactions by R={ij:i,jA,xij>0}. The dynamics of metabolite concentrations S1,,Sn in the network 𝒩 are governed by equations

(10) S˙i=j=1nxjiSj-DiSi,iA

where

(11) Di=j=1nxij,iA.

In what follows, it will be important to distinguish three types of reactions within a module, based on their topological relationship to the I/O metabolites of that module. The topology of the module μ is determined by its set of reactions Rμ={ijR:i,jAμAμIO}. I call all reactions where both the substrate and the product are internal to module μ as reactions internal to μ, and I denote the set of all such reactions by RμiRμ. For example, module μ in Figure 1A has one internal reaction 34. I call all reactions where one of the metabolites is internal to μ and the other is an I/O metabolite as the i/o reactions for μ, and I denote the set of all such reactions by RμioRμ. (I reserve upper-case ‘I/O’ for metabolites and use lower-case ‘i/o’ for reactions.) For example, module μ in Figure 1A has three i/o reactions 13, 14 and 24. Finally, I refer to reactions between any two I/O metabolites for module μ as bypass reactions for module μ, and I denote the set of all such reactions by RμbRμ. For example, module μ in Figure 1A has no bypass reactions but reaction 15 is a bypass reaction for module ν. By definition, all these three sets of reactions Rμi, Rμio and Rμb are non-overlapping, and Rμ=RμiRμioRμb.

Another important concept are the simple paths that lie within a module. For any two metabolites i,jAAμIO, I denote a simple path between them that lies within μ as pijμ or, equivalently as ik1kmj (where all kAμ). I say that each of the metabolites k belongs to (or is contained in) path pijμ (denoted by kpijμ). Similarly, I say that each of the reactions kk+1 belong to (or are contained in) path pijμ (denoted by kk+1pijμ). I will drop superindex μ from pijμ if it is clear what module is being referred to.

Network coarse graining and an algorithm for evaluating the effective rate constant for an arbitrary module

Request a detailed protocol

In this section, I formally introduce and characterize the coarse-graining procedure (CGP). First, I introduce the main idea, which is to eliminate a metabolite that is at QSS and to replace it with a set of new reactions between metabolites adjacent to the eliminated one. This is exactly analogous to the star-mesh transformation in the theory of electric circuits (Versfeld, 1970). The resulting network is a coarse-grained version of the original network in the sense that it has one less metabolite. Next, I define the CGP, which is simply multiple metabolite eliminations applied successively. I prove Proposition 1, which justifies applying the CGP to a whole module and replacing it with a single effective reaction (Corollary 2). Finally, I show how to apply the CGP in practice to compute function F from Equation 1 for modules with some simple topologies.

Elimination of a single metabolite

Request a detailed protocol

I begin by outlining the main idea behind the CGP, which is to replace one metabolite internal to a module with a series of effective reactions between metabolites adjacent to it. If the effective rate constants are defined appropriately, the dynamics of all metabolites in the resulting coarse-grained network are the same as in the original network, provided that the eliminated metabolite is at QSS in the original network.

To formalize this idea, suppose that module μ=(Aμ,xμ) contains m internal metabolites. Let kAμ be the internal metabolites that will be eliminated. Let A{k}=A{k} be the reduced metabolite set and let x{k} be the reduced (n-1)×(n-1) matrix of rate constants defined by 

(12) xij{k}=xij+xikxkjDk,i,jA{k},ij,
(13) xii{k}=0,iA{k},

where Dk is given by Equation 11.

Such metabolite elimination has three properties that follow immediately from Equation 12 and Equation 13. First, the rate constants of reactions do not change as long as the eliminated metabolite does not participate in them. Mathematically, xij{k}=xij for all metabolites i and j that are not adjacent to the eliminated metabolite k. In particular, this is true for all metabolites external to module μ. Second, because equilibrium constants have the property Kij=KiKj for any metabolites i,j,, the rate constants xij{k} obey Haldane’s relationships. Therefore, the reduced metabolite set A{k} and the reduced rate matrix x{k} define a new ‘coarse-grained’ metabolic network 𝒩{k}=(A{k},x{k}). It is easy to show that subnetwork μ after the elimination of metabolite k is still a module. Third, the reaction set of module μ (i.e., its topology) in the coarse-grained network 𝒩{k} depends only on the reaction set of this module in the original network 𝒩, but not on the particular values in the rate matrix xμ.

Next, I will show that the dynamics of metabolites in the coarse-grained network 𝒩{k} are identical to the dynamics of metabolites in the original network 𝒩 where metabolite k is at QSS. Note that if metabolite k is at QSS in the network 𝒩, its concentration is given by

(14) Sk=jAμIOAμ{k}xjkSjDk,

which follows from Equation 10. Now, the dynamics of metabolites in the network 𝒩{k} are governed by equations

(15) S˙i=jA{k}xji{k}SjDi{k}Si,foriA{k},

where Di{k}=jA{k}xij{k}. As mentioned above, xij{k}=xij for all pairs of metabolites where at least one metabolite is external to module μ. Therefore, Equation 15 for the external metabolites are identical to Equation 10 that govern the dynamics of these metabolites in the original network 𝒩. Next, consider the dynamics of the I/O and internal metabolites for module μ in the coarse-grained network 𝒩{k}, that is, those in the set AμIOAμ{k}. For any such metabolite i, the sum in the righthand side of Equation 15 can be re-written as

(16) jAμIOAμ{k}xji{k}Sj=jAμIOAμ{k}(xji+xjkxkiDk)SjxikxkiDkSi=jAμIOAμ{k}xjiSj+xkijAμIOAμ{k}xjkSjDkxikxkiDkSi.

According to Equation 14, the second term in Equation 16 equals xkiSk, so that Equation 16 becomes

(17) jAμIOAμ{k}xji{k}Sj=jAμIOAμxjiSj-xikxkiDkSi.

For any metabolite iAμIOAμ{k}, the second term in the righthand side of Equation 15 can be re-written as

(18) Di{k}=jAμIOAμ{k}(xij+xikxkjDk)-xikxkiDk=Di-xikxkiDk.

Substituting Equation 17 and Equation 18 into Equation 15, we see that Equation 15 is in fact equivalent to Equation 10 for all iA{k} with Sk given by Equation 14.

The coarse-graining procedure (CGP)

Here, I define the CGP for an arbitrary set of internal metabolites by applying metabolite elimination recursively.

Let EAμ be an arbitrary subset of metabolites internal to module μ in the metabolic network 𝒩 and let nE be the number of elements in E. Let AE=AE be the reduced metabolite set after the metabolites have been eliminated. I define the reduced (n-nE)×(n-nE) matrix of rate constants xE as follows. If nE=1, the matrix xE is defined by Equation 12 and Equation 13. If nE>1, then I define it recursively. Suppose that all metabolites in E other than some metabolite kE have been previously eliminated, such that the coarse-grained network 𝒩E=(AE,xE) is defined, with the set of eliminated metabolites E=E{k}, AE=AE, and the known matrix xE. Then, I define the matrix xE through the elimination of metabolite k from 𝒩E, that is, 

(19) xijE=xijE+xikExkjEDkE,i,jAE,ij,
(20) xiiE=0,iAE,

with

(21) DkE=jAExkjE.

I define the the coarse-graining procedure that eliminates the metabolite set E as a map

CGE:𝒩𝒩E=(AE,xE).

As with the elimination of a single metabolite, it is straightforward to show that the rate constants xijE obey Haldane’s relationships, so that 𝒩E is indeed a metabolic network. CGE maps module μ within the metabolic network 𝒩 onto a subnetwork μ within the metabolic network 𝒩E. It is straightforward to show that μ is a module. Whenever there is no ambiguity, I will label both the original and the coarse-grained versions of the module by μ. To simplify notations, if the CGP eliminates the entire module μ (i.e., if E=Aμ), I label it CGμ. I label the coarse-grained network that restults from the application of CGμ by 𝒩μ and I label the effective rate of the reaction substituting module μ in network 𝒩μ as yμ.

Intuitively, the result of coarse-graining should not depend on the order in which the metabolites are eliminated. To see this, let us obtain explicit (i.e. not recursive) expressions for xijE. First, by applying the recursion Equation 19, it is easy to show that elimination of two metabolites E={k,} yields effective rate constants

xij{k,}=xij{k}+xi{k}xj{k}D{k}(22)=xij+Dkxixj+Dxikxkj+xikxkxj+xixkxkjDkDxkxk(23)=xij+xixjD{k}+xikxkjDk{}+xikxkxjDk{}D+xixkxkjD{k}Dk,i,jA{k,},ij.

As expected, Equation 22 and Equation 23 are symmetric with respect to the eliminated metabolites k and . Extrapolating from Equation 23, it is possible to show that for an arbitrary metabolite subset EAμ that contains nE metabolites,

(24) xijE=xij+L=1nE(k1,,kL)xik11xk1k2Dk1E{k1}xkL-1jDkLE{k1,k2,,kL},i,jAE,ij.

Here, the second sum is taken over all nE!/(nE-L)! ordered lists of metabolites (k1,,kL) from E. Each list can be thought of as a simple path within E that connects metabolites i and j. The proof of Equation 24 can be found in Appendix 1. As expected, Equation 24 shows that the effective reation rate xijE does not depend on the order in which metabolites are eliminated. This and other properties of the CGP are listed in Box 1.

Box 1.

Properties of the coarse-graining procedure.

The coarse-graining procedure CGE eliminating metabolites in set EAμ has the following useful properties which follow from Equation 24.

  1. The effective rate constants xijE do not depend on the order in which metabolites are eliminated. Therefore, the composition of coarse-graining procedures is commutative, that is, if E=E1E2, where E1 and E2 are two non-overlapping subsets of Aμ, then

    CGE1CGE2=CGE2CGE1=CGE.
  2. If at least one of the metabolites i or j is not adjacent to any of the eliminated metabolites, then xijE=xij. In particular, xijE=xij if either i and/or j are external to μ.

  3. The topology of module μ after the application of CGE depends only on its original topology but not on the values of its rate constants xμ.

  4. If both metabolites i and j are adjacent to at least one eliminated metabolite, then xijE=xij+α, where α0 depends only on the rate constants of reactions that involve at least one eliminated metabolite. In particular, if both k and are from AE, then xijE is independent of xk.

  5. If E=Aμ, then the effective rate constant yμ of module μ depends on the rate matrix xμ but does not depend on any other reaction rates, that is, Equation 1 holds. Furthermore, the functional form of F(xμ) depends only on the topology of module μ, that is, all modules with the same topology are mapped onto yμ by the same function F.

  6. Suppose that module μ is nested in a larger module ν=(Aν,xν) (see Figure 1A). It follows from Property #1 that CGν=CGμCGAνAμ, that is, yν can be obtained by carrying out the CGP in two stages, by first eliminating module μ and replacing it with the effective reaction with the rate constant yμ and then eliminating the remaining metabolites in Aν. In the network 𝒩μ after the first stage of coarse-graining, all rate constants xνμ are identical to those in the original network, that is, they are independent of xμ, by virtue of Property #2. Therefore yν depends on the rate constants xμ of reactions within module μ only through the effective rate constant yμ of module μ. In other words, Equation 4 holds.

  7. If, in the original network 𝒩, metabolites i and j are adjacent or connected by a simple path that contains only the eliminated metabolites, then metabolites i and j are adjacent in the coarse-grained network 𝒩E.

  8. If, in the metabolic network 𝒩, metabolites i and j are not adjacent (i.e. xij=0) and no simple path exists within the set E (i.e. such that all non-terminal metabolites in this path are from E) that connects metabolites i and j, then metabolites i and j are also not adjacent in the coarse-grained network 𝒩E (i.e., xijE=0).

  9. It follows from properties #7 and #8 that for a simple path p1m=12m to exist in module μ after the application of CGμ, it is neccessary and sufficient that for each i=1,,m-1, either i and i+1 are adjacent in the original module μ or there exists a simple path ii+1 within the original module μ all of whose non-terminal metabolites are from E.

One of key building blocks of the proofs of Theorem 1 and Theorem 2 is the fact that modules can be classified into a finite number of topological classes (see below). To arrive at this classification, it will be convenient to define a composition of coarse-graining procedures, as follows. Suppose that CGE1 and CGE2 are two coarse-graining procedures of network 𝒩 for two subsets of metabolites E1Aμ and E2Aμ. If the sets E1 and E2 are non-overlapping, CGE2 is also defined for the coarse-grained network 𝒩E1 which is the result of applying CGE1 to the original network 𝒩. The result of applying CGE2 to the 𝒩E1 is called the composition of coarse-graining procedures CGE1 and CGE2 of the original network 𝒩 and is denoted as CGE1CGE2.

As defined above, coarse-graining is a formal procedure, and there is no a priori guarantee that (a) it can in fact be carried out for every set of metabolites and (for example, because a metabolite set does not have a steady-state solution); and (b) it will not distort the dynamics in the rest of the network. The following proposition alleviates both of these concerns and thereby justifies the use of the CGP for any subset of internal metabolites within a module (including the entire module μ). It is straightforward to prove it by induction, using the same logic as in the elimination of a single metabolite.

Proposition 1

Request a detailed protocol

Let E be any subset of metabolites internal to module μ. Then,

  1. There exists a joint QSS solution S¯i for all metabolites iE, given the concentrations of the remaining internal and I/O metabolites for module μ.

  2. The dynamics of all remaining metabolites in AE in the coarse-grained metabolic network 𝒩E are the same as in 𝒩 where all metabolites in E are at QSS.

Corollary 1

Request a detailed protocol

Without loss of generality, suppose that the I/O metabolites for module μ are labeled 1 and 2 and its internal metabolites are labeled Aμ={3,4,,m}. There exists a unique QSS S¯i for all iAμ. The QSS concentrations can be obtained by recursively applying equation.

(25) S¯k=1Dk{k+1,,m}(x1k{k+1,,m}S1+x2k{k+1,,m}S2+j{3,,k-1}xjk{k+1,,m}S¯j)

for k=3,4,,m.

Equation 25 follows from Equation 10 for the coarse-grained network obtained by eliminating metabolites k+1,,m.

Corollary 2

Request a detailed protocol

Without loss of generality, suppose that the I/O metabolites for module μ are labeled 1 and 2. Module μ can be replaced with a single effective reaction between its I/O metabolites, whose rate constant yμ can be calculated using Equation 19 and Equation 20 or Equation 24. The dynamics of all metabolites in the resulting coarse-grained metabolic network are identical to their dynamics in the original network N where all metabolites internal to module μ are at the QSS determined by Equation 25.

Computation of effective rate constants for simple modules

Request a detailed protocol

Corollary 2 provides a method for replacing any module μ at QSS with an effective rate yμ=F(xμ), which can be calculated using Equation 19 and Equation 20 or Equation 24. Here, I show how to implement this calculation for three simple metabolic modules.

Linear pathway

Request a detailed protocol

Consider a linear pathway with I/O metabolites 1 and m and internal metabolites 2,3,,m-1 (Figure 6A). This labeling of metabolites is more convenient for the linear pathway. To calculate yμ, I will apply recursion Equation 19 and Equation 20. I start by eliminating metabolite 2. After this initial coarse-graining step, the resulting module is still a linear pathway, where two reactions 123 were replaced with a single reaction 13 with the effective rate constant.

x13{2}=x12x23x21+x23=(1K12x23+1x12)-1.
Simple modules.

(A) Linear pathway. (B) Two parallel pathways.

All other rate constants remain unchanged. Next, I eliminate metabolite 3. The resulting module is still a linear pathway, where now three reactions 1234 were replaced with a single reaction 14 with the effective rate constant

x14{2,3}=x13{2}x34x31{2}+x34=(1K13x34+1K12x23+1x12)-1.

All other rate constants remain unchanged. Continuing this process until all internal metabolites are eliminated, I obtain

(26) yμ=(i=1m-11K1ixii+1)-1,

which is identical to the expression originally obtained by Kacser and Burns, 1973.

Two parallel pathways

Request a detailed protocol

Consider two parallel pathways with I/O metabolites 1 and 2 and internal metabolites 3 and 4 (Figure 6B). I obtain the effective rate constant using Equation 22 with i=1, j=2, k=3, =4. Since x12=x34=0, Equation 22 simplifies to

(27) yμ=D3x14x42+D4x13x32D3D4=x14x42x41+x42+x13x32x31+x32.

Thus, the contributions of parallel pathways are simply added.

Module μ in Figure 1

Request a detailed protocol

To obtain the effective rate constant for module μ shown in Figure 1, I again use Equation 22 with i=1, j=2, k=3, =4.

(28) yμ=D3x14x42+D4x13x32+x13x34x42+x14x43x32D3D4-x34x43,

with D3=x31+x32+x34 and D4=x41+x42+x43.

Classification of modules with respect to ‘marked’ reactions, and the parametric families of functions f1 and f2

The CGP described above allows us to calculate the function F that maps the rate matrix xμ for an arbitrary module μ onto the module’s effective rate constant yμ. F is a multivariate function of the entire matrix xμ. However, in many applications, only one or two reactions are varied at a time while all others remain constant, and we want to know how module’s effective parameter yμ depends on these one or two perturbed reactions. I refer to such singled-out reactions as ‘marked’. When yμ is considered as a function of the rate constant ξ of one marked reaction, I write yμ=f1(ξ), as in Equation 2. When yμ is considered as a function of the rate constants ξ and η of two marked reactions, I write yμ=f2(ξ,η) as in Equation 3.

The functional form of F and, as a consequence, the functional forms of f1 and f2 depend only on the topology of module μ (see Property #5 of the CGP in Box 1). In other words, modules with identical topologies have the same functional forms of f1 and f2, such that each topology of module μ defines a parametric family of functions f1 and f2, where all rate constants within module μ other than ξ, or ξ and η, play a role of parameters.

Since the number of possible topologies is infinite, there is an infinite number of functional forms of F. However, the number of parameteric families of functions f1 and f2 is finite, and it turns out to be small. To see this, notice that for any module with a single marked reaction, the CGP can be carried out in two stages. In the first stage, we can eliminate all metabolites that do not participate in the marked reaction. The resulting coarse-grained module is minimal in the sense that it can have at most two internal metabolites. Such minimal modules (and, as a consequence, all modules with one marked reaction) fall into three distinct topological classes, which are specified by the location of the marked reaction with respect to the I/O metabolites, as shown in Figure 7. This implies that there are only three parameteric families of the function f1. The topologies of the three minimal modules are sufficiently simple that the three corresponding parametric functional forms of f1 can be easily computed by applying the coarse-graining Equation 19 or Equation 22. This result is formulated in Proposition 2.

Classification of single-marked modules.

Left column shows a general module from each topological class. The right column shows a minimal fully connected module in each topological class (see text for details). Circles represent metabolites and lines represent reactions. Only the I/O metabolites and the metabolites that participate in the marked reaction are shown, all other metabolites are suppressed. Short lines that have only one terminal metabolite represent all remaining reactions in which this metabolite participates, reactions between all other metabolites are suppressed. Metabolites are labeled according to the conventions listed in the text. The marked reaction is colored orange and labeled a. The module is represented by a gray rectangle, and the rest of the network is not shown.

The same logic applies to modules with two marked reactions. CGP that eliminates all metabolites that do not participate in the marked reactions maps all such modules onto respective minimal modules, which can have at most four internal metabolites (see Figure 8). This result is formulated in Proposition 3. Minimal modules (and, as a consequence, all modules with two marked reactions) fall into nine distinct topological classes, which are specified by the locations of the marked reactions. All modules from the same topological class are described by functions f2 from the same parametric family. These families are characterized in Corollary 3.

Classification of double-marked modules.

Notations as in Figure 7.

These topological classifications are extremely useful for the following reason. If we can show that all functions from the same parameteric family (corresponding to a given topological class) have some common property irrespectively of the values of the parameters, it would imply that this property holds for all modules from the corresponding topological class. This logic is a key part of the proofs of both Theorem 1 and Theorem 2.

To formalize this reasoning, consider module μ=(Aμ,xμ) and let a=iaja and b=ibjb be two reactions from its set of reactions Rμ. I call a pair (μ,a) a single-marked module and I call a triplet (μ,a,b) a double-marked module, and I refer to reactions a and b as marked within module μ. The topology of a single-marked module (μ,a) is determined not only by the reaction matrix Rμ, but also by the position of the marked reaction, so I refer to the pair (Rμ,a) as the topology of the single-marked module (μ,a). Similarly, I refer to the triplet (Rμ,a,b) as the topology of the double-marked module (μ,a,b). I denote by xμa the matrix of rate constants of all reactions in module μ other than reaction a and I denote by xμ{a,b} the matrix of all rate constants in module μ other than reactions a and b. I denote the sets of all single- and double-marked modules by 1 and 2, respectively. To avoid metabolite labeling ambiguities, I adopt the following conventions:

  1. The I/O metabolites are labeled 1 and 2 and the internal metabolites are labeled 3,4,;

  2. ia,ja{1,2,3,4}; ib,jb{1,2,3,4,5,6};

  3. ia<ja, ib<jb, iaib, jajb.

It is easy to see that the set of all single-marked modules 1 can be partitioned into three non-overlapping topological classes depending on the type of the marked reaction a. I denote the classes of all single-marked modules where the marked reaction is bypass, i/o or internal (see Notations and definitions) by b, io and i, respectively (Figure 7). Similarly, the set 2 can be partitioned into nine non-overlapping topological classes according to the types of marked reactions and the type of metabolite that is shared by both of these reactions (I/O, internal, or none). These classes are listed in Table 1 and illustrated in Figure 8.

Table 1
Classification of double-marked modules.

Metabolites are labeled according to conventions described in the text. m is the minimum number of internal metabolites in a module from class . A is the set of internal and I/O metabolites in all minimal modules in class .

ClassabShared metab.Verbal descriptionmAEquation for f2
b,io,IO(1,2)(1,3)1Bypass and i/o reactions, shared I/O metabolite2{1,2,3}Equation (34)
b,i,(1,2)(3,4)Bypass and internal reactions, no shared metabolies2{1,2,3,4}Equation (35)
io,io,I(1,3)(2,3)3i/o reactions, shared internal metabolite1{1,2,3}Equation (36)
io,io,IO(1,3)(1,4)1i/o reactions, shared I/O metabolite2{1,2,3,4}Equation (37)
io,io,(1,3)(2,4)i/o reactions, no shared metabolites2{1,2,3,4}Equation (38)
io,i,I(1,3)(3,4)3i/o and internal reactions, shared internal metabolite2{1,2,3,4}Equation (39)
io,i,(1,3)(4,5)i/o and internal reactions, no shared metabolites3{1,2,3,4,5}Equation (40)
i,i,I(3,4)(3,5)3Internal reactions, shared internal metabolite3{1,2,3,4,5}Equation (41)
i,i,(3,4)(5,6)Internal reactions, no shared metabolites4{1,2,3,4,5,6}Equation (42)

Each topological class contains infinitely many modules, with various numbers of metabolites and various topologies. However, for each topological class , there is a minimum number of internal metabolites m, such that all modules within must have at least m internal metabolites. I denote the set of metabolites minimal in the topological class by A. It is clear that for the single-marked module classes b and for io, mb=mio=1 and Ab=Aio={3}, and for i, mi=2 and Ai={3,4} (see second column in Figure 7). For the double-marked modules, m and A are given in Table 1 and illustrated in Figure 8.

If a single-marked module from the topological class has the minimum number of metabolites n in that class, I call such module and its topology minimal in . There may be several minimal topologies in a topological class, but there is only one minimal topology that is fully connected. A topology (Rμ,a) is called fully connected if the reaction set Rμ is complete in the sense that it contains reactions between all pairs of metabolites in the minimal metabolite set A. I denote such complete reaction set for the class by R, and I denote the respective fully connected topology by (R,a). I employ the same terminology and analogous notations for double-marked modules. The minimal fully connected topologies are shown in the second column in Figure 7 and Figure 8.

Next, I prove Proposition 2, which is the key step toward the proof of Theorem 1. This proposition states that there are only three functional forms for the function f1 and characterizes them. The idea of the proof is the following. According to Property 5 (Box 1), all single-marked modules that are mapped by the CGP onto a minimal module with the same topology (Rμ,a) have the same functional form of f1. In other words, each minimal topology (Rμ,a) specifies a parameteric family of the function f1. Since the number of possible minimal topologies is finite, the claim of Theorem 1 can be tested for each corresponding functional form of f1. However, the number of minimal topologies is rather large. Fortunately, another simplification is possible. Since the reaction set Rμ of any minimal single-marked module is a subset of the complete reaction set R, the fully connected topology (R,a) specifies the largest parametric family of the function f1 for the class , such that all other families can be obtained from it by setting some parameters to zero, which is equivalent to removing reactions from the fully connected topology. In other words, all single-marked modules that belong to the topological class are described by functions f1 that belong to one parameteric family corresponding to the fully connected topology minimal in . The three parameteric families of f1 are characterized by Proposition 2.

One important consequence of Proposition 2 is that it is not necessary to test the claim of Theorem 1 for each family of f1 that corresponds to each minimal topology. Instead, it is sufficient to test it for the three families that correspond to the fully connected minimal topologies in each class.

Proposition 2

Request a detailed protocol

Let (μ,a) be a single-marked module, and let ξ be the rate constant of reaction a. Then yμ=f1(u), where u=ξ+α for some α0, and the function f1 is given by one of the following expressions.

(29) f1(u)=u,if(μ,a)b,
(30) f1(u)=w12+uw32u/K13+w32,if(μ,a)io,
(31) f1(u)=w12+D3w14w42+D4w13w32+w13w42u+w14w32u/K34D3D4u2/K34, if(μ,a)i.

Here D3=w31+w32+u, D4=w41+w42+u/K34, and all wij0 are independent of ξ.

Proof
Request a detailed protocol

Since any single-marked module (μ,a) belongs to exactly one of three classes b, io and io, I consider these three cases one by one.

Case (μ,a)Mb. According to the labeling conventions outlined above, a=12 (see Figure 7). Equation 29 follows directly from Property #4 of the CGP (Box 1).

Case (μ,a)Mio. According to the labeling conventions, a=13 (see Figure 7). According to Property #1 of the CGP, module μ can be coarse-grained in two stages, by first applying CGAμ{3} which eliminates metabolites 4,,m (those that do not participate in the marked reaction) and then applying CG{3} which eliminates the remaining metabolite 3. Mathematically, CGμ=CGAμ{3}CG{3}. After applying CGAμ{3}, the resulting coarse-grained module μ has a single internal metabolite 3 and at most three effective reactions 12, 13 and 23 (Figure 7), that is, it is minimal in io. By virtue of Properties #2 and #4 of the CGP, the effective rate constants w12, w23 are independent of ξ and uw13=ξ+α. Note that w12 may equal zero, but w230 because μ is a module. Regardless, the reaction set Rμ of module μ is always a subset of the complete reaction set Rio. Thus, to obtain the effective rate constant yμ, I consider the most general case when μ is fully connected and eliminate the remaining internal metabolite 3, which leads to Equation 30.

Case (μ,a)i. According to the labeling conventions, a=34 (see Figure 7). Otherwise, the logic of the proof is exactly the same as for the case (μ,a)io.

Next I prove Proposition 3 which states that, for any double-marked module that belongs to a given topological class, there exists a double-marked module that is minimal in the same class, such that both modules have the same function f2. The corresponding minimal module is obtained from the original module by applying the CGP. This proposition is important because it implies that all functions f2 can be completely characterized by only examing minimal modules. Then, analogously to single-marked modules, Corollary 3 states that function f2 can belong to one of nine parameteric families which are defined by the fully connected minimal topologies in each topological class.

Proposition 3

Request a detailed protocol

Let (μ,a,b) be a double-marked module that belongs to the topological class M, and let ξ and η be the rate constants of reactions a and b, respectively. Then there exist non-negative constants α and β and a double-marked module (μ,a,b) minimal in M such that yμ=yμ=f2(u,v), where

(32) u=ξ+α,
(33) v=η+β

are the rate constants of the marked reactions a and b in μ, respectively, and all other rate constants in μ are independent of ξ or η. Module μ is obtained from μ by the coarse-graining procedure CGμ{a,b} that eliminates all metabolites internal to module μ that do not participate in reactions a or b.

Proof
Request a detailed protocol

To prove this proposition, I will construct the double-marked module (μ,a,b) minimal in by applying CGμ{a,b}. Let m be the mimimal number of internal metabolites in class (see Table 1). According to the metabolite labeling conventions, metabolites n+3,n+4, are neither I/O nor do they participate in the marked reactions. CGμ{a,b} eliminates all these metabolites and maps module μ onto module μ, all of whose internal metabolites participate in reactions a and/or b. Therefore, (μ,a,b) is minimal in class (Figure 8). According to Properties #2 and #4 of the CGP (Box 1), the effective rate constants u and v of reactions a and b in module μ are given by linear relationships in Equation 32 and Equation 33, and the remaining effective rate constants are independent of ξ and η. The fact that yμ=yμ follows from Property #1 of the CGP, CGμ=CGμ{a,b}CGAμ.

Corollary 3

Request a detailed protocol

Let (μ,a,b) be a double-marked module, and let ξ and η be the rate constants of reactions a and b, respectively. The function f2 that maps ξ and η onto module’s effective rate constant yμ belongs to one of nine parametric families. If (μ,a,b)Mb,io,IO, then

(34) f2(u,v)=u+vw32v/K13+w32.

If (μ,a,b)Mb,i,, then

(35) f2(u,v)=u+D3w14w42+D4w13w32+w13w42v+w14w32v/K34D3D4v2/K34,D3=w31+w32+v,D4=w41+w42+v/K34.

If (μ,a,b)Mio,io,I, then

(36) f2(u,v)=w12+uvu/K13+v,

If (μ,a,b)Mio,io,IO, then

(37) f2(u,v)=w12+D3vw42+D4uw32+uw34w42+vw43w32D3D4w34w43,D3=u/K13+w32+w34,D4=v/K14+w42+w43.

If (μ,a,b)Mio,io,, then

(38) f2(u,v)=w12+D3w14v/K24+D4uw32+uw34v/K24+w14w43w32D3D4w34w43,D3=u/K13+w32+w34,D4=w41+v/K24+w43.

If (μ,a,b)Mio,i,I, then

(39) f2(u,v)=w12+D3w14w42+D4uw32+uvw42+w14w32v/K34D3D4v2/K34,D3=u/K13+w32+v,D4=w41+w42+v/K43.

If (μ,a,b)Mio,i,, then

(40) f2(u,v)=W12+W13W32W31+W32,Wij=wij+D4wi5w5j+D5wi4w4j+wi4vw5j+wi5w4jv/K45D4D5v2/K45,D4=w41+w42+w43+v,D5=w51+w52+w53+v/K45,w13u.

If (μ,a,b)Mi,i,I, then

(41) f2(u,v)=W12+D3W14W42+D4W13W32+W13W34W42+W14W43W32D3D4W34W43,Wij=wij+wi5w5jD5,D3=W31+W32+W34,D4=W41+W42+W43,D5=w51+w52+w53+w54,w34u,w35v.

If (μ,a,b)Mi,i,, then

(42)f2(u,v)=W12+D3W14W42+D4W13W32+W13W34W42+W14W43W32D3D4W34W43,(43)Wij=wij+D5wi6w6j+D6wi5w5j+wi5w6jv+wi6w5jv/K56D5D6v2/K56,D3=W31+W32+W34,D4=W41+W42+W43,D5=w51+w52+w53+w54+v,D6=w61+w62+w63+w64+v/K56,w34u.

In Equation 34 through Equation 35, u and v are given by Equation 32 and Equation 33. All effective activities wij0 are constants (other than cases where they stand for u or v) that depend on the topology of module μ and on the parameters xμ{a,b} but do not depend on ξ and η.

Proof
Request a detailed protocol

This statement and Equation 34 through Equation 35 follow directly from Proposition 3 and the fact that the reaction set of any double-marked module in any given topological class is a subset of the complete reaction set in that topological class.

Derivation of Equation 6 and Equation 9

Request a detailed protocol

Consider a higher-level phenotype y, such as the effective activity of a module, which is function of a multivariate lower-level phenotype x=(x1,x2,,xn), such as the rates of individual reactions within the module, y=F(x). Denote the wildtype values of the phenotypes as x0=(x10,x20,,xn0) and y0=F(x0). Consider a mutation that perturbes these values, so that the mutant has lower-level phenotypic values x=(x1,x2,,xn). The relative effect of the mutation on phenotype xi is δxi=xi/xi0-1. If all δxi1 where x denotes the length of vector x, then the value of the higher-level phenotype y in the mutant is given by

(44) y=y0(1+i=1nCiδxi+12i,j=1nHijδxiδxj)+o(δx2).

where

(45) Ci=xi0y0Fxi|x=x0,i=1,,n,
(46) Hij=xi0xj0y02Fxixj|x=x0,i,j=1,,n,

which I refer to as first- and second-order control coefficients of the lower-level phenotypes xi and xj with respect to the higher-level phenotype y.

Now consider two single mutants, A and B, and the double-mutant AB. Each mutation A and B and their combination may perturb all xi phenotypes such that xiA=xi0(1+δAxi), xiB=xi0(1+δBxi), and xiAB=xi0(1+δABxi)=xi0(1+δAxi+δBxi+2δAxiδBxiεxi).

Assuming that δAx1, δBx1 and δABx1, using the approximation in Equation 44 and the definition of εxi (analogous to Equation 5), I obtain

(47) δAy=i=1nCiδAxi+o(δAx),
(48) δBy=i=1nCiδBxi+o(δBx),
(49) εy=i=1nCiεxiδAxiδBxi+12i,j=1nHijδAxiδBxji=1nj=1nCiCjδAxiδBxj+o(1),

where o(1) refers to terms that are vanishingly small as δAxi0 , δBxi0, i=1,n.

I examine two special cases of Equation 49. The first special case is when both mutations affect a single phenotype xk, that is, when all δAxi=0 and all δBxi=0 except for i=k. Then Equation 47, Equation 48, Equation 49 simplify to

(50) δAy=CkδAxk+o(|δAxk|),
(51) δBy=CkδBxk+o(|δBxk|),
(52) εy=εxkCk+Hkk2Ck2+o(1).

Equation 52 is equivalent to Equation 6.

The second special case is when mutation A affects a single phenotypes xk and mutation B affects a single phenotype x (k), i.e., all δAxi=0 except for i=k, all δBxi=0 except for i=, and all εxi=0. Then Equation 47, Equation 48, Equation 49 simplify to

(53) δAy=CkδAxk+o(|δAxk|),
(54) δBy=CδBx+o(|δBx|),
(55) εy=Hk2CkC+o(1).

Equation 55 is equivalent to Equation 9.

Calculation of epistasis in simple modules

Request a detailed protocol

Equation 52 and Equation 55 allow me to compute how epistasis propagates and emerges in arbitrary metabolic networks. In this section, I show how to implement these calculations for three simple metabolic modules considered above and in module ν shown in Figure 3.

Linear pathway

Request a detailed protocol

First, consider how epistasis propagates through a linear pathway (Figure 6A). For simplicity, assume that both mutations A and B affect the same reaction 12. It follows from Equation 26 that

yμ=f1(ξ)=(1ξ+α)-1,

where α is a positive constant. Therefore, the first- and second-order control coefficients of reaction 12 with respect to the flux through the linear pathway μ are given by

C=yμξ=11+αξ,
H=2αξ(1+αξ)2.

Substituting these expressions into the expression for the fixed point ε¯=-H(2C(1-C))-1, I find that ε¯=1, irrespectively of the rates of other reactions in the linear pathway. This implies that epistasis εξ<1 at the level of reaction 12 would induce a lower value of epistasis εyμ<εξ<1 at the level of the entire linear pathway, any value εξ>1 would induce a higher value of epistasis εyμ>εξ>1, and εξ=1 would induce εyμ=1.

Now consider emergence of epistasis in a linear pathway. Suppose that mutation A affects reaction 12 and mutation B affects reaction 23. Denote the rate constant of reactions 12 and 23 by ξx12 and ηx23, respectively. It follows from Equation 26 that

yμ=f2(ξ,η)=(1ξ+1K12η+β)-1,

where β is a positive constant. Therefore,

Cξ=yμξ,Cη=yμK12η,Hξη=2yμ2K12ξη,

which, after substituting into Equation 9, yield εyμ=1. Together with the fact that ε¯=1 (see above), this result implies that epistasis coefficient between any two mutations that affect different reactions in a linear pathway equals 1.

Two parallel pathways

Request a detailed protocol

Suppose that mutation A affects reaction 13 and mutation B affects reaction 14 in the linear metabolic pathway shown in Figure 6B. Denote the rate constants of reaction 13 and 14 by ξx13 and ηx14. It follows from Equation 27 that

yμ=f2(ξ,η)=(1ξ+α)-1+(1η+β)-1,

where α=1/(K13x32) and β=1/(K14x42). Thus, we have Hξη=0, and there is no epistasis between such mutations.

Module ν in Figure 3A

Request a detailed protocol

I denote the rate of the reactions affected by mutations A and B by ξ=x13 and η=x42, and I also denote z=x34. I will calculate the epistasis coefficient εyν in two stages, by first calculating the epistasis coefficient εyμ and then propagating it to εyν using Equation 6. Here I am specifically interested in how εyν depends on the rate constant z.

To compute epistasis between mutations A and B at the level of module μ, I rewrite Equation 28 as

yμ=aξη+bξξ+bηη+cdξη+eξξ+eηη+f,

where a=x14/K13+x32+z, bξ=x32(x41+z/K34), bη=x14(x32+z), c=x14x32z/K34, d=1/K13, eξ=(x41+z/K34)/K13, eη=x32+z, f=x32z/K34+x41z+x32x41. I obtain the following expressions for the first- and second-order control coefficients.

(56) Cξ=ξyμ(c~1z+d~1D)2,
(57) Cη=ηyμ1K14(c~2z+d~2D)2,
(58) Hξη=ξηyμ2zK14(c~1z+d~1)(c~2z+d~2)D3,

where D=dξη+eξξ+eηη+f, c~1=x23/K24+η, d~1=x32(x41+η), c~2=ξ+x14, d~2=x14(ξ/K13+x32). Substituting Equation 56 through Equation 58 into Equation 53 through Equation 55, I obtain

(59) δAyμ=ξyμ(c~1z+d~1D)2δAξ,
(60) δByμ=ηyμ1K14(c~2z+d~2D)2δBη.
(61) εyμ=z(a~z+b~)(c~1z+d~1)(c~2z+d~2),

where a~=c~1c~2 and b~=(ξ/K13+x32)x14η+(x41+η)x32ξ.

To obtain the expression for εyν, I coarse-grain module ν by eliminating the only remaining internal metabolite 2 and obtain

yν=x15+yμx25yμ/K12+x25.

I then apply equation Equation 6 with

(62) C=yμyνx252(yμ/K12+x25)2,
(63) H=2yμ2yνK12x252(u/K12+x25)3.

Figure 3B illustrates how εyν changes as a function of z. It was generated using the following matrix of rate constants:

x=(00.3780.5140.23701.8100001.00142.23200z2.4467.9570z/2.4400.25906.9820.9940.2570).

The Matlab code is available at https://github.com/skryazhi/epistasis_theory.

Next, I consider thress special cases of the toy network depicted in Figure 3A that relate this network to those in Figure 3C and D.

  1. z=0. According to Equation 61, εyμ=0 and hence εyν=H/2C20, with C and H given by Equation 62 and Equation 63. It is easy to see that in this case the reactions 13 and 24 that are affected by mutations are strictly parallel because there is a simple path 1325 that contains only reaction 13 and there is a simple path 1425 that contains only reaction 24 (Figure 3C).

  2. x15=x32=x14=0. In this case, module ν is a linear pathway. Therefore, εyν=1 independently of z, as shown above. This fact can also be obtained directly from Equation 61, Equation 62, Equation 63 and Equation 6.

  3. z. In this case, module μ becomes an effectively linear pathway because most of the metabolic flux between the I/O metabolites 1 and 2 passes through reaction 34. Thus, it follows from Equation 61 that εyμa~/(c~1c~2)=1, as expected. Then, according to Theorem 1, εyν1.

Proof of Theorem 1

As discussed above, the key step toward the proof is Proposition 2, which states that the function f1 belongs to one of three parameteric families, given by Equation 29, Equation 30, Equation 31. To complete the proof, I now explicitly evalute the control coefficient C and the H in Equation 6 for each of these functions and show that the inequalities in Equation 7 and Equation 8 hold for all parameter values.

Proof of Theorem 1

Request a detailed protocol

Let a be the effective reaction within higher-level module ν that represents the lower-level module μ. To simplify notations, I denote yμξ. According to Proposition 2, the functional from of f1 depends only on the topological class of the single-marked module (ν,a). So, I consider the three classes one by one.

Case (ν,a)Mb. From Equation 29, C=ξ/yν and H=0. Therefore, inequalities in Equation 7 and Equation 8 hold.

Case (ν,a)Mio. From Equation 30,

(64) C=ξyν(w32D)2,
(65) H=2ξ2yνw322D31K13=2Cξ/K13D,

where D=(ξ+α)/K13+w32. From Equation 64, it is clear that C0. Re-writing Equation 64 as

C=(ξw32/Dyν)(w32D)

it is also clear that C1 since both ratios in this expression do not exceed 1. From Equation 65 and the fact that 0C1, it follows that ε¯0. To show that ε¯1, note that

D(1-C)=ξ+αK13+w32(1-ξw32/Dyν)ξK13.

Therefore,

ε¯=ξ/K13D(1-C)1.

Therefore, inequalities in Equation 7 and Equation 8 hold.

Case (ν,a)Mi. I re-write Equation 31 as

yν=w12+A~u+B~D,

with u=ξ+α, D=C~u+D~, A~=(w13+w14)(w42+w32/K34), B~=(w31+w32)w14w42+(w41+w42)w13w32, C~=(w31+w32)/K34+(w41+w42), D~=(w31+w32)(w41+w42), which yields

(66) C=ξyνA~D~B~C~D2,
(67) H=2ξ2yν(A~D~B~C~)C~D3=2CC~ξD.

Next, it is straightforward to show that A~D~-B~C~=(w41w32-w31w42)2/K310, which implies that C0. To show that C1, I expand the denominator in Equation 66 and obtain

yνD2(A~u+B~)(C~u+D~)u(A~D~+B~C~)ξ(A~D~-B~C~).

Therefore, numerator in Equation 66 cannot exceed the denominator. The fact that ε¯0 follows directly from Equation 67 together with C1. To show that ε¯1, first note that

yν=w12+A~ξD+A~α+B~DA~ξD.

Therefore,

D(1-C)=C~ξ+C~α+D~(1-A~ξ/Dyν)+ξDB~C~yνC~ξ.

Hence,

ε¯=C~ξD(1C)1.

Therefore, inequalities in Equation 7 and Equation 8 hold in this case as well, which completes the proof.

Proof of Theorem 2

Proving Theorem 2 involves several auxiliary steps. First, I note that any two reactions a and b within any module μ can be either strictly serial, strictly parallel or serial-parallel. Then, Proposition 4 and its Corollary 4 establish that strictly parallel (serial) reactions in (μ,a,b) are also strictly parallel (serial) in a minimal module (μ,a,b), which is obtained from μ by eliminating all metabolites that do not participate in the marked reactions. Next, recall that in both modules (μ,a,b) and (μ,a,b) the same function f2 maps the rate constants of two marked reactions onto module’s effective rate constant (Proposition 3). Since the epistasis coefficient depends only on the shape of this function, we only need to consider all minimal modules in order to understand what kinds of epistasis may arise between mutations affecting strictly serial and strictly parallel reactions in any module. This is a big simplification because the number of different minimal topologies is finite and the parameteric families of function f2 are known for all of them (see Corollary 3). As a consequence, to prove Theorem 2, we could in principle list all of the minimal topologies, identify those where the marked reactions are strictly serial or strictly parallel and evalulate the epistasis coefficient using Equation 9 for every respective function f2.

Unfortunately, the number of minimal topologies is very large, so that such brute-force approach would be quite cumbersome. I take a less cumbersome approach which is based on the realization that a strictly serial or strictly parallel relationship between two reactions cannot be altered by removing a third reaction from the module (Proposition 5). This implies that every minimal topology where the two reactions are strictly serial can be produced from another, more connected, ‘generating’ topology by removing some reactions; and similarly for minimal modules where the reactions are strictly parallel (Proposition 6). All generating topologies can be identified by a simple algorithm given in Appendix 3.

Finally, I prove Theorem 2 in three steps. First, Proposition 7 shows that εyμ0 for any minimal module μ with any strictly parallel generating topology. Second, Proposition 8 shows that that εyμ1 for any minimal module μ with any strictly serial generating topology. Third, the proof of Theorem 2 formally extends this argument to all modules with strictly serial and strictly parallel reactions.

Topological relationships between reactions within a module

Request a detailed protocol

Consider module μ with the I/O metabolites 1 and 2. As described above, a simple path connecting two metabolites i and j within module μ is denoted by pijμ=ikj. If such path contains reactions a1,a2, and does not contain reactions b1,b2,, I denote it as pijμ(a1,a2,,b¯1,b¯2,). I denote the set of all paths pijμ(a1,a2,,b¯1,b¯2,) by 𝒫ijμ(a1,a2,,b¯1,b¯2,).

According to Lemma 1 proven in Appendix 2, any reaction in the module belongs to at least one simple path within module μ that connects the I/O metabolites. Mathematically, 𝒫12μ(a) for any reaction aRμ. Thus, we can define different topological relationships between any two reactions within a module based on the existence of various paths that do or do not contain them. Consequently, we can classify any double-marked module (μ,a,b) based on the toplogocial relationship between its marked reactions. This classification is orthogonal to that given in Table 1.

Two reactions aRμ and bRμ are called parallel within module μ if there exists a simple path p12μ(a,b¯) and a simple path p12μ(b,a¯). Two reactions aRμ and bRμ are called serial within module μ if there exist at least one simple path p12μ(a,b). Two reactions are called strictly parallel within module μ if they are parallel but not serial, they are called strictly serial within module μ if they are serial but not parallel, and they are called serial-parallel within module μ if they are both serial and parallel. It is straightforward to show that there are no other logical possibilities for any two reactions to be anything other than strictly serial, strictly parallel or serial-parallel. This implies that two reactions are strictly parallel if they are not serial, and they are strictly serial if they are not parallel. If reactions a and b are serial, parallel, strictly serial, strictly parallel or serial-parallel within module μ, I also refer to the double-marked module (μ,a,b) as serial, parallel, etc. Since the relationship between reactions depends only on the topology of the module, but not on its rate constants, I also refer to the topology (Rμ,a,b) as serial, parallel, etc.

Recall that coarse-graining procedure CGμ{a,b} eliminates all metabolites internal to module μ other than those participating in reactions a and b. If the double-marked module (μ,a,b) belongs to the topological class , then, according to Proposition 3, CGμ{a,b} maps (μ,a,b) onto a minimal double-marked module (μ,a,b) from the same class . The following proposition, which is easy to prove using Property #9 of the CGP (see Box 1), establishes how this procedure alters the topological relationship between reactions a and b.

Proposition 4

Request a detailed protocol

Let (μ,a,b) be a double-marked module from the topological class M (Table 1) and let (μ,a,b) be the minimal double-marked module in M onto which (μ,a,b) is mapped by CGμ{a,b}.

  1. Reactions a and b are serial in (μ,a,b) if and only if they are serial in (μ,a,b).

  2. If reactions a and b are parallel in (μ,a,b), then they are also parallel in (μ,a,b).

Note that the converse of the second claim in Proposition 4 is not true. In other words, if two reactions a and b are parallel in (μ,a,b), they may not be parallel in (μ,a,b). Figure 9 shows a counter-example illustrating this.

A counter example illustrating that the converse to claim 2 in Proposition 4 may not be true.

Reactions a and b are parallel in (μ,a,b). CGP maps the double-marked module (μ,a,b) onto the minimal double-marked module (μ,a,b) where reactions a and b are not parallel.

Corollary 4

Request a detailed protocol
  1. If reactions a and b are strictly serial in (μ,a,b), they are also strictly serial in (μ,a,b).

  2. If reactions a and b are strictly parallel in (μ,a,b), they are also strictly parallel in (μ,a,b).

  3. If reactions a and b are serial-parallel in (μ,a,b), they are either strictly serial or serial-parallel in (μ,a,b).

Corollary 4 is an important step toward proving Theorem 2. According to Proposition 3, the function that maps the rate constants ξ and η of the reactions a and b in module μ onto the effective rate constant yμ is the same function that maps the rate constants u and v of these reactions in the minimal module μ onto the effective rate constant yμ. It then immediately follows from Equation 9 that the epistasis coefficient εyμ between mutations affecting reactions a and b in the original module μ is the same as the epistasis coefficient εyμ in the minimal module μ. Now, if the reactions a and b are strictly parallel in (μ,a,b), then, according to Corollary 4, these reactions are also strictly parallel in (μ,a,b). Hence, to demonstrate that εyμ0 for any such double-marked module (μ,a,b), it is sufficient to show that εyμ0 for all double-marked modules (μ,a,b) that are minimal in and where the reactions a and b are strictly parallel. And similarly for the strictly serial reactions.

According to this logic, Theorem 2 can be proven by identifying all double-marked modules that are minimal in each of the topological classes listed in Table 1 and where the marked reactions are strictly parallel, evaluating epistasis for all of them, and showing that it is non-positive, irrespectively of the rate constants of any reactions, and similarly for the strictly serial reactions.

Generating topologies

Request a detailed protocol

Since the number of topologically distinct minimal double-marked modules is finite, the approach outlined above is in principle feasible. Unfortunately, the number of topologies to be considered is very large, so in practice it is very cumbersome. To avoid this complication, I take an alternative approach that is based on the same key idea as the proof of Theorem 1. Rather than considering one by one, each minimal topology where the marked reactions are strictly serial or strictly parallel (and the corresponding parametric families of f2), the idea is to identify the most connected minimal topologies (and the corresponding largest parametric families of f2) such that all the other minimal topologies with the strictly serial or strictly parallel reactions (and the corresponding parametric families) can be obtained from them by removing reactions (i.e. setting some parameters to zero).

This idea can be implemented using the following observations. If the two marked reactions are strictly parallel or strictly serial in a minimal module, then removing a third reaction from it does not change this relationship. This statement is proven in Proposition 5. As a consequence, all minimal modules in the topological classes b,io,IO, b,i, and io,io,IO must be strictly parallel because the fully connected minimal topologies are strictly parallel (Figure 8). Similarly, all minimal modules in the topological class io,io,I must be strictly serial because the fully connected minimal topology is strictly serial (Figure 8). The fully connected minimal topologies in all other topological classes are serial-parallel. If the two reactions are serial-parallel, removing a third reaction can change their relationship into a strictly serial or strictly parallel one, depending on which reaction is removed, as shown for example in Figure 3A,C and D. In fact, by removing reactions from the fully connected minimal modules shown in Figure 8, it is easy to show that the topological classes io,io,, io,i,I, io,i,, i,i,I, i,i, contain both strictly serial and strictly parallel modules.

These observations suggests that adding reactions to a minimal module where the marked reactions are strictly serial or strictly parallel will either change their relationship into serial-parallel or will preserve their relationship until the minimal module is fully connected. Therefore, among all minimal modules in a topological class, there must exist the most connected modules where the marked reactions are strictly parallel or strictly serial, such that all other less connected strictly serial or strictly parallel modules can be produced from the most connected ones by removal of reactions. This statement is proven in Proposition 6. Such most connected strictly parallel and strictly serial minimal topologies, which I refer to as ‘generating’, are listed in Table 2 and Table 3. They define the largest parameteric familes of functions f2 which I then examine for the value of εyμ.

Table 2
Strictly parallel generating topologies.
Marked reactionsGenerating topology
ClassabIDExcluded reactionsFigure
b,io,IO1213b,io,IO,FFigure 7
b,i,1234b,i,,FFigure 7
io,io,IO1314io,io,IO,FFigure 7
io,io,1324io,io,,P{34}Figure 9
io,i,I1334io,i,I,P{24}Figure 10
io,i,1345io,i,,P1{34,35}Figure 11
io,i,,P2{25,35}
io,i,,P3{24,25}
i,i,I3435i,i,I,P1{24,25}Figure 12
i,i,I,P2{15,25}
i,i,3456i,i,,P1{35,36,45,46}Figure 13
i,i,,P2{15,16,25,26}
i,i,,P3{24,26,36,45,46}
i,i,,P4{24,25,26,45,46}
i,i,,P5{16,24,25,26,46}
i,i,,P6{14,16,24,26,46}
i,i,,P7{14,15,23,26,35,46}
Table 3
Strictly serial generating topologies.
Marked reactionsGenerating topologies
ClassabIDExcluded reactionsFigure
io,io,I1323io,io,I,FFigure 7
io,io,1324io,io,,S{23}Figure 9
io,i,I1334io,i,I,S1{23}Figure 10
io,i,I,S2{14}
io,i,1345io,i,,S1{14,15}Figure 11
io,i,,S2{23,24,35}
io,i,,S3{15,23,25}
i,i,I3435i,i,I,S1{13,23}Figure 12
i,i,I,S2{23,25,45}
i,i,3456i,i,,S1{23,25,26,45,46}Figure 13
i,i,,S2{13,16,23,26,46}

Proposition 5

Request a detailed protocol

Let (μ,a,b) and (μ,a,b) be two minimal double-marked modules from the same topological class whose sets of reactions are Rμ and Rμ, respectively, and Rμ=Rμ{c} where cRμ{a,b}.

  1. If reactions a and b are strictly parallel in (μ,a,b), they are also strictly parallel in (μ,a,b).

  2. If reactions a and b are strictly serial in (μ,a,b), they are also strictly serial in (μ,a,b).

Proof
Request a detailed protocol

Denote the I/O metabolites in both modules μ and μ by 1 and 2. Since μ and μ are topologically identical except for μ lacking one reaction c, it must be true that 𝒫12μ(a1,a2,,b¯1,b¯2,)𝒫12μ(a1,a2,,b¯1,b¯2,) for any reactions a1,a2,, b1,b2, from Rμ. In other words, there could only be fewer paths connecting the I/O metabolites within module μ compared to module μ. The rest of the proof follows immediately from this fact and the definitions of strictly serial and strictly parallel relatioships.

Next, I define a minimal topology as generating either if it is a fully connected topology (as in topological classes b,io,IO, b,i, and io,io,IO, io,io,I) or if adding any reaction to it would make the marked reactions serial-parallel.

Denote the sets of all double-marked topologies minimal in class where the marked reactions are strictly serial, strictly parallel and serial-parallel by by ser, par and sp, respectively.

Definition 2

Request a detailed protocol

Topology (R,a,b) minimal in M is called a strictly serial generating topology in M if it is strictly serial (i.e. (R,a,b)RMser) and either it is fully connected (i.e. R=RM) or (R{c},a,b)RMsp for any reaction cRMR.

Definition 3

Request a detailed protocol

Topology (R,a,b) minimal in M is called a strictly parallel generating topology in M if it is strictly parallel (i.e. (R,a,b)RMpar) and either it is fully connected (i.e. R=RM) or (R{c},a,b)RMsp for any reaction cRMR.

Clearly, a topological class may have multiple generating topologies, and it is easy to show that every topological class has at least one generating topology. I denote the set of all strictly serial generating topologies for the class by 𝒢ser and I denote the set of all strictly parallel generating topologies for class by 𝒢par. The following proposition justifies the name ‘generating topology’. It states that any strictly serial minimal topology can be produced from some strictly serial generating topology by removing one or multiple reactions, and similarly for any strictly parallel minimal topology.

Proposition 6

Request a detailed protocol

If (R,a,b) is a strictly parallel topology minimal in the topological class M, then there exists a strictly parallel generating topology (Rg,a,b) in M, such that RRg. If (R,a,b) is a strictly serial topology minimal in the topological class M, then there exists a strictly serial generating topology (Rg,a,b) in M, such that RRg.

Proof
Request a detailed protocol

Suppose that is one of the topological classes b,io,IO, b,i,, or io,io,IO. Since the fully connected minimal topology (R,a,b) is strictly parallel, it is a generating topology in . Then, according to Proposition 5, any topology (R,a,b) minimal in is strictly parallel, and Proposition 6 holds. By the same logic, Proposition 6 holds for the topological class io,io,I.

Suppose that is one of the remaining topological classes io,io,, io,i,I, io,i,, i,i,I or i,i,. Then the fully connected minimal topology (R,a,b) is serial-parallel. Suppose that (R,a,b) is strictly parallel. Then R must be a strict subset of R, so that the set C=RR is not empty. Then, either (R,a,b)𝒢par or (R,a,b)𝒢par. If (R,a,b)𝒢par, the Proposition 6 holds. Suppose that (R,a,b)𝒢par. This implies that there exists a reaction c1C, such that R1=R{c1} and (R1,a,b)par ((R1,a,b) cannot be in ser due to Proposition 5). There are three possibilities.

  1. R1=R.

  2. R1R and (R1,a,b)𝒢par.

  3. R1R and (R1,a,b)𝒢par.

Option (a) is not possible since (R1,a,b)par while (R,a,b)sp. Option (b) implies that the Proposition 6 holds. Option (c) implies that there exists a reaction c2C{c1}, such that R2=R1{c2} and (R2,a,b)par, and we have the same three possibilities for R2 as above. Thus, by induction, Proposition 6 must hold. The proof is analogous if (R,a,b) is strictly serial.

Discovering all strictly serial and strictly parallel generating topologies in any given topological class is straightforward because all minimal topologies within can be produced by removing reactions from the unique fully connected topology minimal in shown in Figure 8. In Appendix 3, I provide an algorithm that discovers all strictly serial and strictly parallel generating topologies by sequentially removing reactions from the fully connected minimal topology in each topological class. The code implementing this algorithm in Matlab is available at https://github.com/skryazhi/epistasis_theory. All strictly parallel generating topologies are listed in Table 2 and all strictly serial generating topologies are listed in Table 3. They are also illustrated in Figure 8 and Figure 10 through Figure 14. I label each generating topology by a four-letter combination (see column 4 in Table 2 and Table 3): the first three letters denote the topological class and the last letter (F, P, or S) denotes the specific generating topology within that class. Letter ‘F’ (stands for ‘full’) denotes the fact that the reaction set in the generating topology is complete. Letters ‘P’ (for ‘parallel’) and ‘S’ (stands for ‘serial’) denote strictly parallel and strictly serial generating topologies, respectively; if there are a multiple generating topologies within the same class, they are distinguished by subindices, for example, io,i,,P1; io,i,,P2, etc.

Graphical representation of strictly serial and strictly parallel generating topologies in the class Mio,io,.

Fully connected topology io,io,,F is shown for reference (same as in Figure 8).

Topological relationship between reactions and epistasis

Request a detailed protocol

Each strictly serial and strictly parallel generating topology in a given class (listed in Table 2 and Table 3) is produced by removing reactions from the fully connected topology minimal in (shown in Figure 8). This implies that the parametric family of function f2 that corresponds to any generating topology is obtained from Equation 34 through Equation 35 by setting some parameters wij to zero. In other words, these parametric families are known. Next, I prove Proposition 7 that shows that εyμ0 for every member of every parameteric family of f2 that corresponds to a strictly parallel generating topology and the analogous Proposition 8 for strictly serial topologies.

Now, any minimal strictly parallel topology can in turn be produced by removing reactions from some strictly parallel generating topology, and any minimal strictly serial topology can be produced by removing reactions from some strictly serial generating topology. This implies that any function f2 that corresponds to any strictly parallel minimal topology belongs to the parametric family that corresponds to some strictly parallel generating topology; and any function f2 that corresponds to any strictly serial minimal topology belongs to the parametric family that corresponds to some strictly serial generating topology. Therefore, Propositions 7 and 8 imply that εyμ0 for any minimal strictly parallel topology and that εyμ1 for any minimal strictly serial topology. The proof of Theorem 2 is then concluded by recalling that every strictly parallel module is mapped onto its effective rate constant via function f2 that corresponds to some minimal strictly parallel module, and similarly for strictly serial modules.

Proposition 7

Request a detailed protocol

Let (μ,a,b) be a minimal double-marked module in the topological class M, with u and v being the rates of reactions a and b, respectively, and let y be the effective rate constant of this module. Suppose that mutation A perturbs only reaction a by δAu, and mutation B perturbs only reaction b by δBv, such that |δAu|1, |δBv|1. If reactions a and b are strictly parallel in (μ,a,b), then εy0.

Proof
Request a detailed protocol

According to Proposition 6, the topology of module (μ,a,b) can be produced by removing reactions from some strictly parallel generating topology (Rg,a,b). Therefore, the function f2 that maps u and v in this module onto its effective rate constant y belongs to the parametric family that corresponds to (Rg,a,b). According to Equation 55,

(68) εy=Huv2CuCv+o(1),

where

(69) Cu=uyf2u,
(70) Cv=vyf2v,
(71) Huv=uvy2f2uv.

According to Theorem 1, 0Cu1 and 0Cv1. And since all y>0, u>0, v>0, to prove Proposition 7, it is sufficient to show that 2f2uv0 for any member of any parametric family that corresponds to generating topologies listed in Table 2.

Generating topologies b,io,IO,F and b,i,,F (Figure 8). According to Equation 34 and Equation 35, y=f2(u,v)=u+ϕ(v), which implies that εy=0.

Generating topology io,io,IO,F (Figure 8). According to Equation 37,

y=f2(u,v)=w12+Auv+Bu+BvD,

where D=Euv+Fu+Gv+B, A=w42/K13+w32/K14, B=w32w42+w32w43+w34w42, E=1/(K13K14), F=(w42+w43)/K13, G=(w32+w34)/K14. Therefore,

2f2uv=2w34K14(w32v/K14+B)(w42u/K13+B)D30.

Generating topology io,io,,P (Figure 10). According to Equation 38, y=f2(u,v)=w12+ϕ(u)+ψ(v), which implies εy=0.

Generating topology io,i,I,P (Figure 11). Notice that metabolite 4 together with reactions 14, a and b form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 3 and which is minimal in the topological calss b,io,IO. Denote the effective reaction rate of module μ by y. Therefore, εy=0, as shown above. Since module μ is contained in μ, by Theorem 1, εy0.

Graphical representation of strictly serial and strictly parallel generating topologies in class Mio,i,I.

Fully connected topology io,i,I,F is shown for reference (same as in Figure 8).

Generating topology io,i,,P1 (Figure 12). According to Property 1 of the CGP (Box 1), module μ can be coarse-grained by first eliminating metabolite 3. In the resulting module μ, mutation A perturbs only the rate constant u of the effective reaction a12 (by Properties 2 and 4 of the CGP). Then, according to Equation 50 and Theorem 1, |δAu|1. The double-marked module (μ,a,b) is minimal in the topological class b,i, which implies that εy=0, as shown above.

Graphical representation of strictly serial and strictly parallel generating topologies in class Mio,i,.

Fully connected topology io,i,,F is shown for reference (same as in Figure 8).

Generating topology io,i,,P2 (Figure 12). Module μ can be coarse-grained by first eliminating metabolite 5, which will result in a double-marked module (μ,a,b) that is minimal in the topological class io,io,IO. The rest of the proof for this topology is analogous to that for the topology io,i,,P1.

Generating topology io,i,,P3 (Figure 12). Notice that metabolites 4 and 5 together with reactions a, b, 14, 15, 34 and 35 form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 3 and which is minimal in the topological calss b,i,. The rest of the proof for this topology is analogous to that for the topology io,i,I,P.

Generating topology i,i,I,P1 (Figure 13). Notice that metabolites 4 and 5 together with reactions a, b, 13, 14, 15 and 45 form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 3 and which is minimal in the topological calss io,io,IO. The rest of the proof for this topology is analogous to that for the topology io,i,I,P.

Graphical representation of strictly serial and strictly parallel generating topologies in class Mi,i,I.

Fully connected topology i,i,I,F is shown for reference (same as in Figure 8).

Generating topology i,i,I,P2 (Figure 13). Notice that metabolite 5 together with reactions a, b, and 45 form a double-marked module (μ,a,b) whose I/O metabolites are 3 and 4 and which is minimal in the topological calss b,io,IO. The rest of the proof for this topology is analogous to that for the topology io,i,I,P.

Generating topology i,i,,P1 (Figure 14). According to Equation 35, y=f2(u,v)=x12+ϕ(u)+ψ(v), which implies εy=0.

Graphical representation of strictly serial and strictly parallel generating topologies in class Mi,i,.

Fully connected topology i,i,,F is shown for reference (same as in Figure 8).

Generating topology i,i,,P2 (Figure 14). Notice that metabolites 5 and 6 together with reactions a, b, 35, 36, 45 and 56 form a double-marked module (μ,a,b) whose I/O metabolites are 3 and 4 and which is minimal in the topological calss b,i,. The rest of the proof for this topology is analogous to that for the topology io,i,I,P.

Generating topology i,i,,P3 (Figure 14). Module μ can be coarse-grained by first eliminating metabolites 4 and 6, which will result in a double-marked module (μ,a,b) that is minimal in the topological class io,io,IO. The rest of the proof for this topology is analogous to that for the topology io,i,,P1.

Generating topology i,i,,P4 (Figure 14). Module μ can be coarse-grained by first eliminating metabolite 4, which will result in a double-marked module (μ,a,b) that is minimal in the topological class io,i, with a strictly parallel generating topology io,i,,P3. The rest of the proof for this topology is analogous to that for the topology io,i,,P1.

Generating topology i,i,,P5 (Figure 14). Module μ can be coarse-grained by first eliminating metabolite 6, which will result in a double-marked module (μ,a,b) that is minimal in the topological class i,i,I with a strictly parallel generating topology i,i,I,P1. The rest of the proof for this topology is analogous to that for the topology io,i,,P1.

Generating topology i,i,,P6 (Figure 14). Notice that metabolites 4 and 6 together with reactions a, b, 35, 36, 45 form a double-marked module (μ,a,b) whose I/O metabolites are 3 and 5 and which is minimal in the topological calss io,io,. The rest of the proof for this topology is analogous to that for the topology io,i,I,P.

Generating topology i,i,,P7 (Figure 14). Using Equation 42, I show in Appendix 4 that

(72) 2f2uv=2βK31(Auv+Bu)(Avu+Bv)(Euv+Fu)3,

where Av, Bv, Eu, Fu are all non-negative constants and β0.

Proposition 8

Request a detailed protocol

Let (μ,a,b) be a minimal double-marked module in the topological class M, with u and v being the rates of reactions a and b, respectively, and let y be the effective rate constant of this module. Suppose that mutation A perturbs only reaction a by δAu, and mutation B perturbs only reaction b by δBv, such that |δAu|1, |δBv|1. If reactions a and b are strictly serial in (μ,a,b), then εy1.

Proof
Request a detailed protocol

The logic of the proof is the same as for Proposition 7, that is, it is sufficient to show that εy1 for any double-marked module (μ,a,b) with any strictly serial generating topology listed in Table 3.

Generating topology io,io,I,F (Figure 8). According to Equation 36,

(73) y=f2(u,v)=w12+uvD

where D=u/K13+v. Therefore,

Cu=(vD)2uy,Cv=1K12(uD)2vy,Huv=2K121yD(uvD)2.

Substituting these expressions into Equation 68, I obtain

εy=yuv/D1

because yuv/D according to Equation 73.

Generating topology io,io,,S (Figure 10). According to Property 1 of the CGP (Box 1), module μ can be coarse-grained by first eliminating metabolite 3. In the resulting module μ, mutation A perturbs only the rate constant u of the effective reaction a14 (by Properties 2 and 4 of the CGP). Then, according to Equation 50 and Theorem 1, |δAu|1. The double-marked module (μ,a,b) is minimal in the topological class io,io,I which implies that εy1, as shown above.

Generating topology io,i,I,S1 (Figure 11). Notice that metabolite 3 together with reactions a, b, and 14 form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 4 and which is minimal in the topological calss io,io,I. Therefore, if the effective reaction rate of module μ is y, εy1, as shown above. According to Equation 50, Equation 51 and Theorem 1, |δAy|1, |δBy|1. Since module μ is contained in μ, by Theorem 1, εy1.

Generating topology io,i,I,S2 (Figure 11). Module μ can be coarse-grained by first eliminating metabolite 4, which results in a double-marked module (μ,a,b) that is minimal in the topological class io,io,I. The rest of the proof for this topology is analogous to that for the topology io,io,,S.

Generating topology io,i,,S1 (Figure 12). Module μ can be coarse-grained by first eliminating metabolites 4 and 5, which results in a double-marked module (μ,a,b) that is minimal in the topological class io,io,I. The rest of the proof for this topology is analogous to that for the topology io,io,,S.

Generating topology io,i,,S2 (Figure 12). Notice that metabolites 3 and 4 together with reactions a, b, 14, 15 and 34 form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 5 and which is minimal in the topological calss io,io, with the strictly serial generating topology io,io,,S. The rest of the proof for this topology is analogous to that for the topology io,i,I,S1.

Generating topology io,i,,S3 (Figure 12). Notice that metabolites 3 and 5 together with reactions a, b, 14, 34, and 35 form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 4 and which is minimal in the topological calss io,io, with the strictly serial generating topology io,io,,S. The rest of the proof for this topology is analogous to that for the topology io,i,I,S1.

Generating topology i,i,I,S1 (Figure 13). Notice that metabolite 3 together with reactions a, b, and 45 form a double-marked module (μ,a,b) whose I/O metabolites are 4 and 5 and which is minimal in the topological calss io,io,I. The rest of the proof for this topology is analogous to that for the topology io,i,I,S1.

Generating topology i,i,I,S2 (Figure 13). Notice that metabolites 3 and 5 together with reactions a, b, 13, 14, and 15 form a double-marked module (μ,a,b) whose I/O metabolites are 1 and 4 and which is minimal in the topological calss io,i,I with the strictly serial generating topology io,i,I,S2. The rest of the proof for this topology is analogous to that for the topology io,i,I,S1.

Generating topology i,i,,S1 (Figure 14). Module μ can be coarse-grained by first eliminating metabolites 5 and 6, which results in a double-marked module (μ,a,b) that is minimal in the topological class io,i,I with the strictly serial generating topology io,i,I,S1. The rest of the proof for this topology is analogous to that for the topology io,io,,S.

Generating topology i,i,,S2 (Figure 14). Module μ can be coarse-grained by first eliminating metabolite 6, which results in a double-marked module (μ,a,b) that is minimal in the topological class i,i,I with the strictly serial generating topology i,i,I,S1. The rest of the proof for this topology is analogous to that for the topology io,io,,S.

Proof of Theorem 2

Request a detailed protocol

According to Proposition 3, the coarse-graining procedure CGμ{a,b} maps the double-marked module (μ,a,b) onto a double-marked module (μ,a,b) that is minimal in the same topological class as (μ,a,b), and the rates u, v of reactions a, b in μ are given by linear relations in Equation 32 and Equation 33. Clearly, |δAu|1 and |δBv|1. Furthermore, none of the other reaction rates wij in μ depend on ξ or η, so that δAwij=0 and δBwij=0 for all wij other than u and v, and εwij=0 for all wij including u and v. It then follows from Proposition 3 that εyμ=εyμ.

Now, according to Corollary 4, if reactions a and b are strictly parallel in (μ,a,b), they are also strictly parallel in (μ,a,b). Therefore, by Proposition 7, εyμ0. Analogously, if reactions a and b are strictly serial in (μ,a,b), they are also strictly serial in (μ,a,b). Therefore, by Proposition 8, εyμ1.

Sensitivity of Theorem 1 and Theorem 2 with respect to the magnitude of mutational effects

Request a detailed protocol

According to Proposition 2, function f1 for any module belongs to one of three parametric families, which correspond to the three minimal fully connected modules shown in Figure 7. As mentioned in the Results, for modules in the class b, function f1 is linear, so that the claims of Theorem 1 continue to hold for mutations with finite effects. To evaluate the sensitivity of Theorem 1 with respect to the effect sizes of mutations for the topological classes io and i, I generated 1000 minimal single-marked modules ν from each of these topological classes with random parameters. Evaluating only minimal modules is sufficient because for any module from a given topological class there exists a minimal module from the same class, such that both of them map the lower level phenotype yμ onto the higher-level phenotype yν via the same function f1 (see Proposition 2).

To this end, I drew each xij (i<j) from a mixture of a point measure at 0 (with weight 0.25) and an exponential distribution with mean 1 (with weight 0.75). The point measure at 0 ensures that minimal modules that are not fully connected are represented in the sample. I drew each Kij (i<j) as a ratio of two random numbers from an exponential distribution with mean 1. As a result, the distribution of non-zero xij values had the interdecile range of (5.7×10-2,3.91) with median 0.65.

I denote the effective rate constant of the reaction that represents the lower-level module μ by ξyμ. In modules from the topological class io, it is reaction 13 and in modules from the topological class i, it is the reaction 34. I perturbed ξ by two mutations A and B with relative effects δAξ and δBξ and epistasis εξ. I chose nine different pairs of mutational effects (δAξ,δBξ):(0.01,0.01), (-0.1,-0.1), (-0.5,-0.5), (0.01,0.01), (0.1,0.1), (0.5,0.5), (-0.01,0.01), (-0.1,0.1), (-0.5,0.5), and 16 different values of εξ ranging from −1 to 2 with an increment of 0.2. Since the rate constant ξAB of the double mutant cannot be negative, I skipped those combinations of perturbations and epistasis values for which δAξ+δBξ+2(εξ)(δAξ)(δBξ)<-1. I then computed the resulting values δAyν, δByν and εyν at the level of the effective rate constant yν of the higher-level module ν.

Using these data, I inferred the function ϕ that maps lower-level epistasis εξ onto higher-level epistasis εyν, as follows. For any minimal single-marked module from the topological classes io or i, the effective rate constant yν can be written as

yν=x12+A~ξ+B~D,

where D=C~ξ+D~ and A~=x32, B~=0, C~=1/K13 , D~=x32 for modules from the topological class io (see Equation 30), and A~=(x13+x14)(x42+x32/K34), B~=(x31+x32)x14x42+(x41+x42)x13w32, C~=(x31+x32)/K34+(x41+x42) , D~=(x31+x32)(x41+x42) for modules from the topological class i (see Equation 31). Therefore, for any perturbation δξ, we have

δyν=A~D~-B~C~D2ξyνδξ1+(C~ξ/D)δξ.

Since δABξ is a linear function of εξ, δAByν is a hyperbolic function of εξ. Therefore, εyν is also a hyperbolic function of εξ,

(74) εyν=ϕ(εξ)=a-bεξ+c,

where constants a, b and c depend on the parameters of module ν and on the mutational effect sizes δAξ and δBξ. I numerically calculated these parameters for each sampled module and each pair of mutational effects.

The main results of Theorem 1 are that, when the effects of mutations are infinitesimal, the map ϕ has a fixed point ε¯, this fixed point is located between 0 and 1, and it is unstable. I use equation Equation 74 to test whether these statements also hold when the effects of mutations are finite. Specifically, it is easy to see that the map ϕ has a fixed point ε¯ if the discriminant d=(a-c)2-4(b-ac) is positive. In this case, I designate ε¯ as the one of two roots 1/2(a-c±d) that is closer to zero. I then check whether this fixed point is located between 0 and 1. I check whether it is unstable by comparing the derivative of ϕ at ε¯ with 1.

According to Proposition 6, function f2 for any module where the reactions affected by mutations are strictly parallel belongs to one of 17 parameteric families, which correspond to the strictly parallel generating topologies listed in Table 2. And similarly, function f2 for any module where the reactions affected by mutations are strictly serial belongs to one of 11 parameteric families, which correspond to the strictly serial generating topologies listed in Table 3. Therefore, to evaluate the sensitivity of Theorem 2 with respect to the effect sizes of mutations I generated 104 double-marked modules (μ,a,b) with each of the strictly serial and strictly parallel topologies with random parameters. I drew xij and Kij as described above. I chose the same nine pairs of mutational effects (δAξ,δBη) as above, where ξ and η are the rate constants of reactions affected by mutations A and B: (-0.01,-0.01), (-0.1,-0.1), (-0.5,-0.5), (0.01,0.01), (0.1,0.1), (0.5,0.5), (-0.01,0.01), (-0.1,0.1), (-0.5,0.5).

I found that, for some modules, individual mutational perturbations δAyμ and/or δByμ at the level of the whole module were too small, which resulted in numerical instabilities. To avoid them, I calculated epistasis εyμ only for cases where the effects of both mutations δAyμ and δByμ exceeded the precision threshold of 10-5. As a result, I evaluated epistasis in less than 104 modules per generating topology and pair of mutational effects, but this number never fell below 1000. When comparing the values of epistasis with 0 and 1, I used the same precision threshold of 10-5 to avoid numerical problems. In addition, I found that for mutations affecting strictly serial reactions there is a substantial fraction of modules where εyμ falls between 0.99 and 1 (see Figure 4—figure supplement 3). This is not a numerical artifact, but probably reflects real clustering of epistasis coefficients around 1, which is expected for the linear pathway irrespective of its parameters (see above).

The Matlab code for this analysis is available at https://github.com/skryazhi/epistasis_theory.

Kinetic model of glycolysis

I downloaded the kinetic metabolic model of E. coli glycolysis by Chassagnole et al., 2002 from the BioModels database (Malik-Sheriff et al., 2019) on September 15, 2015 (model ID BIOMD0000000051). I used the Matlab SimBiology toolbox to interpret the model. To validate the model, I simulated it for 40 s and reproduced Figures 4 and 5 from Chassagnole et al., 2002. The Matlab code is available at https://github.com/skryazhi/epistasis_theory.

Modifications to the original model

Request a detailed protocol

I simplified and modified the model by (a) fixing the concentrations of ATP, ADP, AMP, NADPH, NADP, NADH, NAD at their steady-state values given in Table V of Chassagnole et al., 2002 and (b) removing dilution by growth. I then created four models of sub-modules of glycolysis by retaining the subsets of metabolites and enzymes shown in Figure 5—figure supplement 1 and Table 4 and removing other metabolites and enzymes. Each sub-module has one input and one output metabolite. Note that, since some reactions are irreversible, it is important to distinguish the input metabolite from the output metabolite. The concentrations of the input and the output metabolites in each model are held constant at their steady-state values given in Table 4. I defined the flux through the sub-module as the flux toward the output metabolite contributed by the sub-module (Table 4). This flux is the equivalent of the quantitative phenotype yμ of a module in the analytical model. In addition, I made the following modifications specific to individual sub-modules.

Table 4
Definition of modules in the glycolysis network shown in Figure 5—figure supplement 1.

Enzyme abbreviations are listed in Table 6. Metabolite abbreviations are listed in Table 5.

ModelInternal metabolitesConcentrations of I/O metabolitesReactionsOutput flux
UGPP6 pg, dhap, e4p, f6p, fdp, rib5p, ribu5p, sed7p, xyl5p[g6p]=3.82 mM, [gap]=0.44 mMALDO, G6PDH, PFK, PGDH, PGI, Ru5P, R5PI, TA, TIS, TKa, TKbJALDO+JTIS+JTKb+JTKa-JTA
LG2 pg, 3 pg, pgp[gap]=0.44 mM, [pep]=0.08 mMENO, GAPDH, PGK, PGMJENO
GPPall in UGPP and in LG, gap[g6p]=3.82 mM, [pep]=0.08 mMall in UGPP and in LGJENO
FULLall in GPP, g6p, pep[Ext glu]=2 µM, [pyr]=10 µMall in GPP, PTS, PK, PEPCxylJPK+JPTS
Table 5
Names of metabolites used in the kinetic model of glycolysis.
2 pg2-Phosphoglycerate
3 pg3-Phosphoglycerate
6 pg6-Phosphogluconate
dhapDihydroxyacetonephosphate
e4pErythrose-4-phosphate
f6pFructose-6-phosphate
fdpFructose-1,6-bisphosphate
g6pGlucose-6-phosphate
gapGlyceraldehyde-3-phosphate
gluGlucose
pepPhosphoenolpyruvate
pgp1,3-Diphosphoglycerate
pyrPyruvate
rib5pRibose-5-phosphate
ribu5pRibulose-5-phosphate
sed7pSedoheptulose-7-phosphate
xyl5pXylulose-5-phosphate
Table 6
Names of enzymes used in the kinetic model of glycolysis.
ALDOAldolase
ENOEnolase
G6PDHGlucose-6-phosphate dehydrogenase
GAPDHGlyceraldehyde-3-phosphate dehydrogenase
PFKPhosphofructokinase
PGDH6-Phosphogluconate dehydrogenase
PGIGlucose-6-phosphateisomerase
PGKPhosphoglycerate kinase
PGMPhosphoglycerate mutase
PEPCxylPEP carboxylase
PKPyruvate kinase
PTSPhosphotransferase system
R5PIRibose-phosphateisomerase
Ru5PRibulose-phosphate epimerase
TATransaldolase
TISTriosephosphate isomerase
TKaTransketolase, reaction a
TKbTransketolase, reaction b
  1. In the FULL model, the stoichiometry of the PTS reaction was changed to

    [Ext glu]+[pep][g6p]+[pyr]

    and the value of the constant KPTS,a1 was set to 0.02 mM, based on the values found in the literature (Stock et al., 1982; Natarajan and Srienc, 1999).

  2. In all models other than FULL, the extracellular compartment was deleted.

  3. In all models, the concentrations of the I/O metabolites were set to values shown in Table 4, which are the steady-state concentrations achieved in the FULL model with the concentration of extracellular glucose being 2 µM and pyruvate concentration being 10 µM.

Calculation of flux control coefficients and epistasis coefficients

Request a detailed protocol

I calculate the first- and second-order flux control coefficients (FCC) Ci and Hij for flux J with respect to reactions i and j as follows (see Equation 45 and Equation 46). I perturb the rmax,i of reaction i by factor between 0.75 and 1.25 (10 values in a uniformly-spaced grid), such that δrmax,i[-0.25,0.25]. Then, I obtain the steady-state flux J in each perturbed model and calculate the flux perturbations δJ=J/J0-1, where J0 is the corresponding flux in the unperturbed model. Then, to obtain Ci and Hii, I fit the linear model

δJCi(δrmax,i)+Hii2(δrmax,i)2

by least squares. If the estimated value of Ci was below 10-4 for a given sub-module, I set Ci to zero and exclude this reaction from further consideration in that sub-module because it does not affect flux to the degree that is accurately measurable. If the estimated value of Hii is below 10-4, I set Hii to zero.

To calculate the non-diagonal second-order control coefficients Hij, I create a 4×4 grid of perturbations of δrmax,i and δrmax,j and calculate the resulting flux perturbations δJ (16 perturbations total). Since Ci, Cj, Hii and Hjj are known, I obtain Hij, by regressing

δJ-(Ci(δrmax,i)+Hii2(δrmax,i)2)-(Cj(δrmax,j)+Hjj2(δrmax,j)2)

against

(δrmax,i)(δrmax,j).

If the estimated value of Hij is below 10-4, I set Hij to zero. I estimate the epistasis coefficient εJ between mutations affecting reactions i and j as

εJ=Hij2CiCj.

Establishing the topological relationships between pairs of reactions

Request a detailed protocol

To establish the topological relationship (strictly serial, strictly parallel, or serial-parallel) between two reactions, I consider the smallest module (LG, UGPP, GPPP, or FULL) which contains both reactions. I then manually identify whether there exists a simple path connecting the input metabolite with the output metabolite for that module that passes through both reactions. (Note that, since some reactions are irreversible in this model, it is important to distinguish the input metabolite from the output metabolite). If such path does not exist, I classify the topological relationship between the two reactions as strictly parallel. If such path exists, I check if there are two paths connecting the input to the output metabolites such that each path contains only one of the two focal reactions. If such paths do not exist, I classify the topological relationship between the two reactions as strictly serial. Otherwise, I classify it as serial-parallel.

Appendix 1

Proof of Equation 24

First, the terms in Equation 24 can be re-arranged as follows 

(75) xijE=xij+xik1xk1jDk1E{k1}+xik2xk2jDk2E{k2}++xik1xk1k2xk2jDk1E{k1}Dk2E{k1,k2}++xik2xk2k1xk1jDk2E{k2}Dk1E{k1,k2}+=xij+xik1Dk1E{k1}(xk1j+xk1k2xk2jDk2E{k1,k2}+)=xk1jE{k1}+xik2Dk2E{k2}(xk2j+xk2k1xk1jDk1E{k1,k2}+)=xk2jE{k2}+=xij+ExixjE{}DE{}.

Next, I will demonstrate the validity of Equation 75 by induction. It is clear that for E={k}Aμ, Equation 24 reduces to Equation 12. Now suppose that nE>1 and that there exists a metabolite kE, such that Equation 75 holds for the subset E=E{k}, that is, 

(76) xijE=xij+ExixjE{}DE{}.

I will now show that then Equation 75 also holds for E. To do so, I use the definition of xijE (Equation 19) and apply Equation 76 to expand terms xijE and xikE as follows.

(77) xijE=xijE+xikExkjEDkE=xij+ExixjE{}DE{}+xikxkjEDkE+ExixkE{}xkjEDE{}DkE=xij+ExiDE{}(xjE{}+xkE{}xkjEDkE)+xikxkjEDkE.

Recalling that E=E{k}, I re-write Equation 75 as 

(78) xijE=xij+ExixjE{}DE{}+xikxkjEDkE.

Since the first and third terms of Equation 77 and Equation 78 are identical, to complete the proof, it is sufficient to show that for any E,

(79) 1DE{}(xjE{}+xkE{}xkjEDkE)1DE{k,}(xjE{k,}+xkE{k,}xkjE{k}DkE{k})=xjE{}DE{}.

To show that Equation 79 holds, I first use Equation 19 and Equation 21 to express xkjE{k} and DkE{k} in terms of effective reaction rates after the elimination of the metabolite set E{k,}

(80) xkjE{k}=xkjE{k,}+xkE{k,}xjE{k,}DE{k,}DkE{k}=DkE{k,}xkE{k,}xkE{k,}DE{k,},

which imply that

1DE{k,}(xjE{k,}+xkE{k,}xkjE{k}DkE{k})=DkE{k,}DE{k,}DkE{k}(xjE{k,}+xkE{k,}xkjE{k,}DkE{k,})=DkE{k,}xjE{}DE{k,}DkE{k}.

Finally, it follows from Equation 80 that DE{k,}DkE{k}=DkE{k,}DE{}, which completes the proof of Equation 79. Thus, Equation 75 and equivalently Equation 24 hold for any metabolite subset EAμ.

Appendix 2

Existence of a simple path that contains a given reaction

Lemma 1

Let μ be a module with the reaction set Rμ. Then, for any reaction aRμ, there exists a simple path p12(a) within μ that connects the I/O metabolites and contains reaction a, that is, P12μ(a).

Proof

Reaction a is either a bypass, i/o, or internal reaction for module μ. If a is a bypass reaction, then the statement is trivially true. If a is an i/o reaction, then, without loss of generality, let a=1j. Since μ is a module, there exists a simple path jj12 that connects the internal metabolite j to the I/O metabolite 2. Therefore, the path 1jj12 connects the I/O metabolites and contains reaction a.

Suppose that a=ij is an internal reaction. To prove the statement, it is sufficient to show that there exists a pair of non-intersecting paths p1i and p2i, such that one of them contains a and the other does not. Since μ is a module, there exists a pair of non-intersecting paths p1i and p2i and a pair of non-intersecting paths p1j and p2j within module μ (I omitted super-index μ to simplify notations). There are two mutually exclusive possibilities. (i) Metabolite j is contained in either of the paths p1i or p2i and/or metabolite i is contained in either of the paths p1j or p2j. (ii) Metabolite j is not contained in either of the paths p1i or p2i and metabolite i is not contained in either of the paths p1j or p2j, that is, jpui and ipuj, u=1,2. It is trivial to construct the neccessary paths p1i and p2i in case (i).

Appendix 2—figure 1
Illustration for the proof of Lemma 1.

Consider case (ii). If paths p2j and p1i do not intersect, then let p1i=p1i and

p2i=ij=a2,=p2j

and the statement is true. Suppose paths p2j and p1i intersect. Then, among all metabolites that belong to both p1i and p2j, let metabolite k be the one closest to j along the path p2j (see Figure). Then the segment pkj of path p2j and the path p1i do not intersect. Let 

p1i=1kalong p1i along p2jji=a,

If p1i′′ and p2i do not intersect, then the lemma is true. If p1i′′ and p2i do intersect, this intersetion can only occur within the segment pkj of path p1i′′, excluding metabolites k and j (see Figure). This is because the remaining segment p1k of path p1i′′ is also a segment of p1i, which, by assumption, does not intersect p2i. Suppose that among all metabolites that belong to both the segment pkj of path p1i′′ and the path p2i metabolite is the one closest to j along the path p1i′′. Then let

p1i=p1i,p2i=2along p2ialong p2jji=a.

The path p2i does not intersect the path p1i because its first segment p2 belongs to path p2i and its second segment pj belongs to the segment pkj of path p1i′′ (and, as mentioned above, segment pkj does not intersect p1i). Thus, the statement holds for case (ii) as well.

Appendix 3

An algorithm for discovering all strictly serial and strictly parallel generating topologies

Suppose that (R,a,b)𝒢ser, that is, (R,a,b) is a strictly serial generating topology in the topological class . Since RR where R is the complete reaction set R for class , R can be discovered by sequentially removing reactions from R. The same logic holds for strictly parallel generating topologies. The following algorithm implements this idea.

  1. Define function generate_topology_list. This function takes a topology (R,a,b)sp as input and returns a new list of topologies L as output, which is produced as follows. Initialize L=. For every reaction ciR{a,b}, construct the reaction subset Ri=R{ci} and use Definition 1 to test whether Ri corresponds to a valid module. If Ri corresponds to a module, add (Ri,a,b) to list L; otherwiese, discard. It can be proven that, as long as (R,a,b)sp, there exists at least one ciR, such that Ri corresponds to a module, that is, L. Return list L.

  2. Initialization.

    1. Pick a topological class .

    2. Test whether (R,a,b)ser. If so, 𝒢ser={(R,a,b)} and 𝒢par=. Return 𝒢ser, 𝒢par.

    3. Test whether (R,a,b)par. If so, 𝒢par={(R,a,b)} and 𝒢ser=. Return 𝒢ser, 𝒢par.

    4. Set 𝒢ser=, 𝒢par=. Use function generate_topology_list with (R,a,b) as input and obtain the list of reaction sets L. Proceed to Step 3 with list L.

  3. Take list L=((R1,a,b),(R2,a,b),,(Rk,a,b)) as input. Set L=. Proceed to Step a with i=1.

    1. Test whether (Ri,a,b) belongs to par, ser or sp. Choose one of the alternatives b, c, or d.

    2. (Ri,a,b)par. If (Ri,a,b) is not in the set 𝒢par and if (Ri{c},a,b)sp for all cRRi, then add (Ri,a,b) to 𝒢par. Proceed to Step e.

    3. (Ri,a,b)ser. If (Ri,a,b) is not in the set 𝒢ser and if (Ri{c},a,b)sp for all cRRi, then add (Ri,a,b) to 𝒢ser. Proceed to Step e.

    4. (Ri,a,b)sp. Use function generate_topology_list with (Ri,a,b) as input and obtain the list of reaction sets Li. Replace L with LLi. Proceed to Step e.

    5. If i=k, proceed to Step f. Otherwise, proceed to Step a with i+1.

    6. If L, then proceed to Step 3 with L as input. Otherwise, return 𝒢ser, 𝒢par.

Appendix 4

Derivation of Equation 72

In this case, Equation (42) simplify to

W12(v)=w12+w16w52v/K56D56(v),W13(v)=w13+w16w63D5(v)D56(v),W14(v)=w16w54v/K56D56(v),W23(v)=w25w63vD56(v),W24(v)=w24+w25w54D6(v)D56(v),W34(u,v)u+w36w54v/K56D56(v),

where

D56(v)=D5(v)D6(v)v2K56,D5(v)=w52+w54+v,D6(v)D6(v)=w61+w63+v/K56.

Notice that the only effective rate constant that depends on u is W34, and W34u=1. Thus, it is easy to differentiate y, given by Equation 35, with respect to u, if we isolate the term W34 in both the numerator and the denomintor, 

(81) y=W12+A56,34W34+B56,34D56,34,

where

D56,34=E56,34W34+F56,34,A56,34=W14W42+W13W32K34+W13W42+W14W32K34(W13+W14)(W32K34+W42),B56,34=W14W42(W31+W32)+W13W32(W41+W42),E56,34=W31+W32K34+(W41+W42),F56,34=(W31+W32)(W41+W42).

It is also useful to obtain another expression for y, which is easier to differentiate with respect to v. To do that, we can first eliminate metabolites 3 and 4 to obtain effective reaction rates

V12(u)=w12+w13w42uD34(u),V15(u)=w13w45uD34(u),V16(u)=w16+w13w36D4(u)D34(u),V25(u)=w25+w24w45D3(u)D34(u),V26(u)=w24w36u/K34D34(u),V56(u,v)=v+w36w54u/K34D34(u),

where

D34(u)=D3(u)D4(u)u2K34,D3(u)=w31+u+w36,D4(u)=w42+u/K34+w45.

The only effective activity that depends on v is V56, and V56v=1. Thus, isolating the term V56 (and recalling that V65=V56/K56), we obtain the following expression for y, which is easy to differentiate with respect to v.

(82) y=V12+A34,56V65+B34,56D34,56,

where

D34,56=E34,56V65+F34,56,A34,56=(V16+V15)(V62K65+V52),B34,56=V15V52(V61+V62)+V16V62(V51+V52),E34,56=V61+V62K65+(V51+V52),F34,56=(V61+V62)(V51+V52).

Using symbolic computation it is possible to show that (see Mathematica notebook [Supplementary file 1] and Mathematica notebook pdf [Supplementary file 2], p. 2)

(83) D56,34D56=D34,56D34.

Also notice that, any module with the reaction set i,i,,P7 is symmetric with respect to swapping metabolite labels 1 with 2, 3 with 5, and 4 with 6. It is easy to check that Equations (81), (82) respect this symmetry.

Differentiating Equation 81 with respect to u, after some algebra I obtain

(84) yu=yW34=1K31(W31W42-W32W41D56,34)2.

Analogously, differentiating Equation 82 with respect to v, I obtain

(85) yv=1K56yV65=1K51(V52V61-V51V62D34,56)2.

Notice that Equation 85 can also be obtained from Equation 84 by symmetry with respect to the aforementioned metabolite relabeling.

Next, using symbolic computation (see Mathematica notebook (Supplementary file 1) and Mathematica notebook pdf (Supplementary file 2), p. 3), it is possible to show that

(86) (W31W42-W32W41)D56D56,34D56=Auv+BuEuv+Fu,

where all coefficients

Au=w31K56ψ+w42ϕ,Bu=ψϕ,Eu=u[w31w52+w31w54+w36w52K34K56+(w42w61+w42w63+w45w61)]+D4(u)ϕ+D3(u)K56ψ,Fu=u(w52+w54K34ϕ+(w61+w63)ψ)+ϕψ.

and

ϕ=w31w61+w31w63+w36w61,ψ=w42w52+w42w54+w45w52

are independent of u and v and are non-negative. Similarly (see Mathematica notebook (Supplementary file 1) and Mathematica notebook pdf (Supplementary file 2), p. 4),

(87) (V52V61-V51V62)D34D34,56D34=Avu+BvEvu+Fv,

where

Av=w61ψ+w52K34ϕ,Bv=ψϕ,Ev=v[(w42w61+w42w63+w45w61)+w31w52+w31w54+w36w52K34K56]+D5(v)K34ϕ+D6(v)ψ,Fv=v((w42+w45)ϕ+w31+w36K56ψ)+ϕψ.

We can now obtain the second derivative 2yuv, taking into account Equation 86. Alternatively, we can obtain 2yuv by differentiating yv with respect to u, taking into account Equation 87. The denominators in both expressions would be identical due to Equation 83. Therefore the expression for the second derivative must have the form given by Equation 72, that is,

2yuv=2βK31(Auv+Bu)(Avu+Bv)(Euv+Fu)3,

where β is independent of u and v. Thus, according to Equation 84, Equation 86,

β=AuFu-BuEuAvu+Bv=-(w36ψ/K56+w45ϕ),

which is verified in Mathematica notebook (Supplementary file 1) and Mathematica notebook pdf (Supplementary file 2), p. 4. Similarly, according to Equation 85, Equation 87,

βK31=1K51AvFvBvEvAuv+Bu=w63ψ+w54/K34ϕK51

which is verified in Mathematica notebook (Supplementary file 1) and Mathematica notebook pdf (Supplementary file 2), p. 5.

Appendix 5

Relationship to the flux balance analysis

Here, I discuss the relationship of the model presented in this paper to the flux balance analysis (FBA), a widely used approach to modeling whole-cell metabolism (Orth et al., 2010). FBA and my model are designed to address different questions. My model was designed to explore how the flux depends on the kinetic parameters of a module. FBA was designed to avoid this dependence. Nevertheless, the two models are conceptually and mathematically related and in fact are in some sense equivalent, as discussed below. The most important similarity is that both models rely on flux balance—the assumption that internal metabolites are at steady state—as a key simplification.

The major difference is that my model assumes that all reaction kinetics are first order, which allows for analytical tractability, while FBA is agnostic with respect to reaction kinetics and takes into account only their stoichiometries and possibly some additional ‘capacity’ constraints on the reaction rates (Orth et al., 2010). This difference deserves some additional discussion. In real cells, the concentrations of internal metabolites and the resulting fluxes (which are functions of these concentrations) are determined by the stoichiometries of reactions, the activities of enzymes and the concentrations of external nutrients. A general approach to modeling such systems, adopted by the metabolic control analysis and related theories, is to explicitly specify the dynamic equations for the metabolites, such as Equation 10, and solve them to find steady-state concentrations of internal metabolites (see e.g. Savageau, 1976). The steady-state fluxes are then determined automatically. In my model, the internal steady state exists and is unique for any module (see Corollary 1). In contrast, there are no dynamic equations within FBA. Instead, the steady-state fluxes are subject only to the mass-balance equations and the capacity conditions, which typically form a severely underdetermined system. Therefore, arriving at a unique flux distribution in FBA requires additional assumptions. The typical approach is to first fix at least some nutrient uptake rates and then maximize an objective function, such as the growth rate or biomass yield (Feist and Palsson, 2010). In other words, FBA trades-off the ability to handle reactions with arbitrary kinetics and unknown kinetic parameters against the necessity to impose auxiliary conditions to find the right solution among many plausible ones.

As a consequence of these model design choices, mutations that add or remove reactions can be naturally studied within the FBA framework, but there is no natural way to incorporate mutations that perturb the kinetic parameters of reactions (He et al., 2010; Alzoubi et al., 2019). In contrast, mutations of small or large effect, including reaction additions and deletions, can be in principle naturally studied within my model.

Due to differences in the assumptions, my model and FBA describe different sets of biological systems. However, there are special cases which can be described by both FBA and my model, and in these special cases the two models are mathematically equivalent.

From the perspective of FBA, modules with two I/O metabolites and first-order kinetics described by my model are a special case. So, let us apply FBA to a metabolic module μ=(Aμ,xμ) with two I/O metabolites 1 and 2. For the purposes of FBA, let the I/O metabolite 1 be the ‘external nutrient’ and let the I/O metabolite 2 be the ‘biomass’. Suppose that metabolites in the set AinAμ are adjacent to the external nutrient and metabolites in the set AoutAμ are adjacent to the biomass. Then reactions 1i for iAin, that is, those that convert the nutrient into intermediate metabolites, are the ‘uptake reactions’. And the reactions i2 for iAout, that is, those that convert intermediate metabolites into biomass, are the ‘biomass reactions’. Denote the rate of reaction consuming metabolite i and producing metabolite j by vij and let Jin=iAinv1i and Jout=iAoutvi2 be the input and output fluxes, respectively. To study epistasis within the FBA framework, we would need to obtain Jout.

The assumption of my model that all reaction kinetics are first-order translates into the FBA formalism as the fact that the elements of the stoichiometry matrix S take values −1, +1 or 0 and there are no capacity constraints on the rates of reactions. Then, the mass-balance equation Sv=0 leads to the equality Jin=Jout, so that the flux through the module does not depend on module’s internal structure. Instead, it must be given as an auxiliary condition.

In my model, the steady-state flux J through the module depends in general on the internal structure of the module (through the effective rate constant yμ) and is given by J=yμ(S1-S2/K12) for any concentrations S1 and S2 of the external nutrient and the biomass. However, in the degenerate case where all uptake reactions are irreversible (i.e., K1i= for all iAin), J=Jin=iAinx1iS1, such that it is independent of the internal structure of the module. In this sense, this special case of my model is equivalent to FBA. However, my model and FBA are not equivalent in terms of the distribution of internal fluxes. My model still produces a unique steady sate and the corresponding flux distribution, whereas FBA does not specify a unique flux distribution in this case without additional auxiliary conditions.

Appendix 6

Connection between metabolism and growth

Here I describe a simple ‘bioreactor’ model for connecting a hierarchical metabolic network described in this paper with cellular growth. Suppose that module μ with I/O metabolites 1 and n is at the top level of the metabolic hierarchy and describes the whole cell. The I/O metabolite 1 can be thought of as nutrients, and the I/O metabolite n can be thought of as proteins, so that cellular metabolism converts nutrients into biomass. Denote the concentrations of metabolites and proteins within cells by Si and p. Suppose that cells have a fixed volume v and the number of cells is N. Then, the absolute abundance of proteins and metabolites across all N cells is P=pvN and Ai=SivN, i=1,,n-1.

The proteins produced by the cell are the enzymes that catalyze all the reactions inside module μ. Let the relative expression level of the enzyme catalyzing reaction ij, i<j be ϕij (with i=1n-1j>iϕij=1), so that its concentration inside cells is ϕijp. If the specific forward and reverse activities of this enzyme are aij and aji (so that aij obey the Haldane equalities and aii=0), then the total forward and reverse activities are xij=aijϕijp and xji=ajiϕijp.

Finally, I assume that all reactions converting internal metabolites into biomass are irreversible, and I assume that the cell density in the bioreactor is small enough that the nutrient concentration S1 stays constant. In other words, I model early exponential growth. Then the dynamics of metabolite and protein abundances are governed by equations

(88) A˙i=i=1n1xjiAjDiAi,i=2,n1.
(89) P˙=i=2n1xinAi,

where Di=i=1nxij.

I assume that all cells are identical and at steady state, such that protein and metabolite concetrations inside cells are constants and the production of proteins and metabolites manifests itself in the multiplication of cells. Then Equation 88 and Equation 89 become

(90) λSi=i=1n1xjiSjDiSi,i=2,n1.
(91) N˙=λN,

where λ=i=2n-1xinSi is the steady-state flux into biomass which defines the exponential growth rate of the system.

Finally, assuming that the rates of consumption and production of all metabolites are much greater than their dilution by growth, I set λSi0. Then Equation 90 defines the same steady state for module μ as described in the section Network coarse-graining. Therefore, module μ can be replaced with the effective activity xμ between nutrients and biomass, yielding λ=xμ.

Data availability

All data generated or analyzed during this study are included in the manuscript and supporting files. Code is available on GitHub.

References

  1. Book
    1. Fisher RA
    (1930)
    The Genetical Theory of Natural Selection
    Oxford: The Clarendon Press.
  2. Conference
    1. Kacser H
    2. Burns JA
    (1973)
    The control of flux
    Symposia of the Society for Experimental Biology.
    1. Kacser H
    2. Burns JA
    (1981)
    The molecular basis of dominance
    Genetics 97:639–666.
    1. Keightley PD
    (1989)
    Models of quantitative variation of flux in metabolic pathways
    Genetics 121:869–876.