Emergence and propagation of epistasis in metabolic networks
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 wholeorganism phenotypes, such as fitness. Here, I derive two rules for how epistasis between mutations with small effects propagates from lower to higherlevel phenotypes in a hierarchical metabolic network with firstorder 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 intergene epistasis should be common, and it should generically depend on the genetic background and environment. Furthermore, the epistasis coefficients measured for highlevel 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 higherorder 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 largescale 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 higherorder 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 ‘highlevel’ phenotype, such as fitness, which depends on the performance of the focal ‘lowerlevel’ 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 lowerlevel phenotypes to higherlevel 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, BagheriChaichian 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, longterm 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 highlevel phenotypes, such as fitness. Segrè et al., 2005 and He et al., 2010 used flux balance analysis (FBA, Orth et al., 2010) to compute genomewide 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 genomescale 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 lowerlevel phenotypes onto a higherlevel 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 firstorder 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 lowerlevel phenotypes manifests itself at higherlevel wholeorganism phenotypes, such as fitness, and what kind of information may be gained from measurements of such higherlevel 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,\mathrm{\dots},n\}$ with concentrations ${S}_{1},\mathrm{\dots},{S}_{n}$ which can be interconverted by reversible firstorder biochemical reactions. The rate of the reaction converting metabolite $i$ into metabolite $j$ is ${x}_{ij}\left({S}_{i}{S}_{j}/{K}_{ij}\right)$ where ${K}_{ij}$ is the equilibrium constant. The rate constants ${x}_{ij}$, which satisfy the Haldane relationships ${x}_{ji}={x}_{ij}/{K}_{ij}$ (CornishBowden, 2013), form the matrix $\overrightarrow{x}={\parallel {x}_{ij}\parallel}_{i,j=1}^{n}$. The metabolite set $A$ and the rate matrix $\overrightarrow{x}$ define a biochemical network $\mathcal{N}=(A,\overrightarrow{x})$.
The firstorder 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 ${x}_{ij}$ depend on the concentrations and the specific activities of enzymes and therefore can be altered by mutations. ${K}_{ij}$ characterize the fundamental chemical nature of metabolites $i$ and $j$ and cannot be altered by mutations (Savageau, 1976).
The wholecell 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 graphtheoretic sense) if there exists an enzyme that catalyzes a biochemical reaction between them, that is, if ${x}_{ij}>0$. Now consider a subset of metabolites ${A}_{\mu}\subset A$. For this subset, let ${A}_{\mu}^{\mathrm{IO}}$ be the set of all metabolites that do not belong to ${A}_{\mu}$ but are adjacent to at least one metabolite from ${A}_{\mu}$. Let ${\overrightarrow{x}}_{\mu}$ be the submatrix of $\overrightarrow{x}$ which corresponds to all reactions where both the product and the substrate belong to either ${A}_{\mu}$ or ${A}_{\mu}^{\mathrm{IO}}$. The metabolite subset ${A}_{\mu}$ and the rate matrix ${\overrightarrow{x}}_{\mu}$ form a subnetwork $\mu =({A}_{\mu},{\overrightarrow{x}}_{\mu})$ of network $\mathcal{N}$. I refer to ${A}_{\mu}$ and ${A}_{\mu}^{\mathrm{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.
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. nonselfintersecting) 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 $i\mathrm{\in}{A}_{\mu}$, 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; BagheriChaichian 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}_{\mu}={y}_{\mu}\left({S}_{1}{S}_{2}/{K}_{12}\right)$, where
is the effective reaction rate constant of module μ (Figure 1B). Importantly, ${y}_{\mu}$ depends only on the rate matrix ${\overrightarrow{x}}_{\mu}$, 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 $\mathcal{N}$ can be coarsegrained by replacing module μ at steady state with a single firstorder biochemical reaction with rate ${y}_{\mu}$. Importantly, such coarsegraining 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 starmesh 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}_{\mu}$ quantifies the function of module μ (a macroscopic parameter) while the rates ${\overrightarrow{x}}_{\mu}$ describe the specific biochemical implementation of the module (microscopic parameters).
The effective rate constant ${y}_{\mu}$ of module μ depends on the entire rate matrix ${\overrightarrow{x}}_{\mu}$. 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}_{\mu}$ depends on the one or two rate constants of the perturbed reactions. When ${y}_{\mu}$ is considered as a function of the rate constant ξ of one reaction, I write
and when ${y}_{\mu}$ is considered as a function of the rate constants ξ and η of two reactions, I write
The rate constants of all other reactions within module μ play a role of parameters in functions $f}_{1$ and $f}_{2$.
Consider now a network $\mathcal{N}$ that has a hierarchical structure, such that there is a series of nested modules $\mu \subset \nu \subset \mathrm{\cdots}$, in the sense that ${A}_{\mu}\subset {A}_{\nu}\subset \mathrm{\cdots}$ (Figure 1A). Since any module at steady state can be replaced with an effective firstorder biochemical reaction, there exists a hierarchy of quantitative metabolic phenotypes ${y}_{\mu},{y}_{\nu},\mathrm{\dots}$ (Figure 1B,C). These phenotypes are of course functionally related to each other. Specifically, because ν is a ‘higherlevel’ module (in the sense that it contains a ‘lowerlevel’ module μ), the matrix ${\overrightarrow{x}}_{\nu}$ can be decomposed into two submatrices ${\overrightarrow{x}}_{\mu}$ and ${\overrightarrow{x}}_{\nu \setminus \mu}$ where the latter is the matrix of rate constants of reactions that belong to module ν but not to module μ. Since replacing the lowerlevel module μ with an effective reaction with rate constant ${y}_{\mu}$ does not alter the dynamics of metabolites outside of μ, ${y}_{\nu}$ must depend on all elements of ${\overrightarrow{x}}_{\mu}$ only through ${y}_{\mu}$, that is,
where rates ${\overrightarrow{x}}_{\nu \setminus \mu}$ act as parameters of function $f}_{1$. Thus, in the hierarchy of metabolic phenotypes ${y}_{\mu},{y}_{\nu},\mathrm{\dots}$, a phenotype at each subsequent level depends on the phenotype at the preceding level according to Equation 4, and the lowest level phenotype ${y}_{\mu}$ 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 ${x}_{ij}$, such that the wildtype value ${x}_{ij}^{0}$ changes to ${x}_{ij}^{A}$. This mutation can be quantified at the microscopic level by its relative effect ${\delta}^{A}{x}_{ij}={x}_{ij}^{A}/{x}_{ij}^{0}1$. If the reaction between metabolites $i$ and $j$ belongs to nested modules $\mu ,\nu ,\mathrm{\dots}$, then mutation A may impact the functions of these modules, which can be quantified by the relative effects ${\delta}^{A}{y}_{\mu}$, ${\delta}^{A}{y}_{\nu}$, etc. at each level of the hierarchy.
Consider now another mutation B that only perturbs the rate constant ${x}_{k\mathrm{\ell}}$ 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 higherlevel 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}_{\mu}$ 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
where ${\delta}^{AB}{y}_{\mu}$ denotes the effect of the combination of mutations A and B on phenotype ${y}_{\mu}$ relative to the wildtype. Since I only consider two mutations A and B, I will write $\epsilon {y}_{\mu}$ instead of ${\epsilon}^{AB}{y}_{\mu}$ to simplify notations. Note that other epistasis coefficients can always be computed from $\epsilon {y}_{\mu}$, ${\delta}^{A}{y}_{\mu}$ and ${\delta}^{B}{y}_{\mu}$, 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 $\epsilon {y}_{\mu}$ at a lower level μ of the metabolic hierarchy, what can we say about their epistasis coefficient $\epsilon {y}_{\nu}$ 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 ${x}_{ij}$ of one enzyme and mutation B only perturbs the activity ${x}_{k\mathrm{\ell}}$ of another enzyme that belongs to the same module μ, then what values of $\epsilon {y}_{\mu}$ 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 lowerlevel μ are small, it follows from Equation 4 and Equation 5 that
where $C={f}_{1}^{\prime}{y}_{\mu}/{y}_{\nu}$ and $H={f}_{1}^{\prime \prime}{y}_{\mu}^{2}/{y}_{\nu}$ are the first and secondorder control coefficients of the lowerlevel module μ with respect to the flux through the higherlevel module ν and $o\left(1\right)$ 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 $\varphi $ with slope $1/C$ and a fixed point $\overline{\epsilon}=H{\left(2C\left(1C\right)\right)}^{1}$, which both depend on the topology of the higherlevel module ν and the rate constants ${\overrightarrow{x}}_{\nu \setminus \mu}$.
To gain some intuition for how the map $\varphi $ 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 $f}_{1$ is monotonically increasing (i.e. the higher ${y}_{\mu}$, the more flux can pass through the linear pathway ν) and concave (i.e. as ${y}_{\mu}$ grows, other reactions in ν become increasingly more limiting, such that further gains in ${y}_{\mu}$ yield smaller gains in ${y}_{\nu}$). Indeed, it is easy to show that $C={\left(1+\alpha {y}_{\mu}\right)}^{1}>0$ and $H=2\alpha {y}_{\mu}{\left(1+\alpha {y}_{\mu}\right)}^{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 $\epsilon {y}_{\mu}$ that already arose at the lower level would propagate to negative epistasis $\epsilon {y}_{\nu}$ 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 $\epsilon {y}_{\mu}$ was already sufficiently large at the lower level, it would induce even larger positive epistasis $\epsilon {y}_{\nu}$ at the level of the linear pathway ν. In fact, when module ν is a linear pathway, $\overline{\epsilon}=1$, so that $\epsilon {y}_{\nu}>1$ whenever $\epsilon {y}_{\mu}>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).
Theorem 1
For any module ν,
and
The proof of Theorem 1 is given in Materials and methods. Its main idea is the following. The functional form of $f}_{1$ 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 $f}_{1$. However, this is not the case. In fact, all higherlevel modules that contain a lowerlevel module fall into three topological classes defined by the location of the lowerlevel module with respect to the I/O metabolites of the higherlevel module (see Proposition 2 and Figure 7 in Materials and methods). To each topological class corresponds a parametric family of the function $f}_{1$, 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 $\varphi $ from epistasis at a lowerlevel to epistasis at the higherlevel 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, $\epsilon {y}_{\nu}\le \epsilon {y}_{\mu}<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, $\epsilon {y}_{\mu}>1$, it will induce even stronger positive epistasis at the next level of the hierarchy, that is, $\epsilon {y}_{\nu}\ge \epsilon {y}_{\mu}>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<\epsilon {y}_{\mu}<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 $\epsilon {y}_{\mu}$, the topology of the higherlevel module ν and the microscopic rate constants ${\overrightarrow{x}}_{\nu \setminus \mu}$.
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 lowerlevel phenotype cannot turn into positive epistasis for a higherlevel 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 $\epsilon {y}_{\nu}$ can take values in all three domains described above, depending on the biochemical details of this module. Using the recursive procedure for evaluating ${y}_{\mu}$ described in Materials and methods, it is straightforward to obtain an analytical expression for ${y}_{\nu}$ as a function of the rate matrix ${\overrightarrow{x}}_{\nu}$, from which $\epsilon {y}_{\mu}$ can also be obtained (see Materials and methods for details). To demonstrate that $\epsilon {y}_{\mu}$ 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 $z\equiv {x}_{34}$ 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 $\epsilon {y}_{\mu}$ varies as a function of $z$ for one particular choice of all other rate constants. When $z$ is small, $\epsilon {y}_{\mu}<0$. As $z$ increases, it becomes weakly positive ($0<\epsilon {y}_{\mu}<1$) and eventually strongly positive ($\epsilon {y}_{\mu}>1$). Thus, in my model, there are no fundamental constraints on the types of epistasis that can emerge between mutations.
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, ‘higherorder’ and ‘environmental’ epistasis are generic features of metabolic networks.
Upon closer examination, the toy example in Figure 3 also suggests that the sign of $\epsilon {y}_{\nu}$ 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. nonselfintersecting) 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 serialparallel. 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 serialparallel.
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 ${\delta}^{A}\mathit{}\xi $ and mutation B perturbs only the other reaction by ${\delta}^{B}\mathit{}\eta $. In the limit ${\delta}^{A}\mathit{}\xi \mathrm{\to}\mathrm{0}$ and ${\delta}^{B}\mathit{}\eta \mathrm{\to}\mathrm{0}$, the following statements are true. If the affected reactions are strictly parallel then $\epsilon \mathit{}{y}_{\mu}\mathrm{\le}\mathrm{0}$. If the affected reactions are strictly serial, then $\epsilon \mathit{}{y}_{\mu}\mathrm{\ge}\mathrm{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
where ${C}_{\xi}=\frac{\partial {f}_{2}}{\partial \xi}\frac{\xi}{{y}_{\mu}}$, ${C}_{\eta}=\frac{\partial {f}_{2}}{\partial \eta}\frac{\eta}{{y}_{\mu}}$, ${H}_{\xi \eta}=\frac{{\partial}^{2}{f}_{2}}{\partial \xi \partial \eta}\frac{\xi \eta}{{y}_{\mu}}$ are the first and secondorder control coefficients of the affected reactions with respect to the flux through module μ and $o(1)$ denotes terms that vanish when ${\delta}^{A}\xi $ and ${\delta}^{B}\eta $ 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 $\epsilon {y}_{\mu}$ for an arbitrary module μ, we need to know the first and second derivatives of function $f}_{2$. Analogous to function $f}_{1$, there is a finite number of parametric families to which $f}_{2$ 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 $f}_{2$ (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 serialparallel. And it is easy to show that not all members of each topological class have the same sign of $\epsilon {y}_{\mu}$. However, modules from the same topological class where the affected reactions are strictly parallel or strictly serial fall into a finite number of topological subclasses (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 subfamilies of function $f}_{2$. For all members of these subfamilies, Equation 9 yields $\epsilon {y}_{\mu}\le 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 subfamilies of function $f}_{2$. For all members of these subfamilies, Equation 9 yields $\epsilon {y}_{\mu}\ge 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 higherlevel modules that contain a lowerlevel module fall into three topological classes, which I label ${\mathcal{M}}^{\mathrm{b}}$, ${\mathcal{M}}^{\mathrm{io}}$ and ${\mathcal{M}}^{\mathrm{i}}$, depending on the location of the lowerlevel module within the higher level module (see Figure 7). The topological class specifies the parametric family of the function $f}_{1$ which maps the effective rate constant ${y}_{\mu}$ onto the effective rate constant ${y}_{\nu}$ (see Equation 4). For all modules from the topological class ${\mathcal{M}}^{\mathrm{b}}$, function $f}_{1$ 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 ${\mathcal{M}}^{\mathrm{io}}$ and ${\mathcal{M}}^{\mathrm{i}}$, function $f}_{1$ 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 nonlinear function $\varphi $ that maps the epistasis coefficient $\epsilon {y}_{\mu}$ onto the epistasis coefficient $\epsilon {y}_{\nu}$ 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 lowerlevel module ${y}_{\mu}$ 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}_{\mu}$, map $\varphi $ had the same properties as described in Theorem 1, even for mutations with large effect, that is, it had a fixed point $\overline{\epsilon}$ in the interval $[0,1]$ and this fixed point was unstable. When the effects of both mutations on ${y}_{\mu}$ 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 $\overline{\epsilon}$ shifted slightly above 1. As the magnitude of mutational effects increased, the fraction of sampled modules with $\overline{\epsilon}>1$ grew, reaching 42% when both mutations increased ${y}_{\mu}$ by 50%. In most of these cases, $\overline{\epsilon}$ remained below 2, and I found only one module with $\overline{\epsilon}>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 $\varphi $ was very close to the identity map. For 289 modules (14.5%), the fixed point disappeared when both mutations increased ${y}_{\mu}$ by 50%. In all these cases, $\epsilon {y}_{\nu}<\epsilon {y}_{\mu}$, indicating that even large positive epistasis may decline as it propagates through the metabolic hierarchy when the effects of mutations are finite.
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 $f}_{2$ which maps the rate constants ξ and η of the affected reactions onto the effective rate constant ${y}_{\mu}$. To test how well Theorem 2 holds when the effects of mutations are finite, I calculated $\epsilon {y}_{\mu}$ 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 higherlevel 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 firstorder 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 secondorder, reaction kinetics are nonlinear, and in several cases the reaction rates are modulated by other metabolites (Chassagnole et al., 2002).
Testing the predictions of the analytical theory in this computational model faces two complications. First, in a nonlinear 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 nonlinearly 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}_{\nu}$ 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 steadystate 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 nonzero 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 nonoverlappng 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 lowerlevel phenotype cannot turn into positive epistasis for a higherlevel 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 glucose6phosphate isomerase (PGI) and 6phosphogluconate 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 higherlevel 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 highlevel 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 (BagheriChaichian 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 firstorder kinetics. Networks with only firstorder 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 firstorder 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 coarsegraining 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 highlevel 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 genomewide 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 lowerlevel phenotype they also generally exhibit epistasis for a higherlevel phenotype. The second contribution comes from the fact that lowerlevel phenotypes map onto higherlevel phenotypes via nonlinear 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 lowerlevel phenotype, epistasis must emerge for the higherlevel 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 downregulated 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\times G\times G$ interactions (higheroder epistasis) and $G\times G\times 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 intergene 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 fitnessrelated phenotypes among genomewide 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 highlevel 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 higherlevel 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 nonexclusive 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 highlevel 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 coarsegrained 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 withinprotein 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 nonlinear 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 nonpleiotropic.
If epistasis is in fact generic for highlevel phenotypes, why do we not observe it more frequently? For example, a recent study that tested almost all pairs of gene knockout 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 intergene epistasis are interpreted. In particular, contrary to a common belief, a nonzero 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 protocolBefore 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 quasisteady 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 ‘coarsegrain’ the network by removing metabolites. This approach is a standard biochemicalnetwork reduction technique (Segel, 1988); for example, the BriggsHaldane derivation of the MichaelisMenten formula is based on this idea. In general, the resulting effective reactions have more complex (nonlinear) kinetics than the original reactions. However, when the original reactions are firstorder, the effective reactions are also firstorder, 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 coarsegraining 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}_{\mu}$ of this reaction. In other words, the CGP is an algorithm for evaluating function $F({\overrightarrow{x}}_{\mu})$ 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}_{\nu}$ of the larger module is to first coarsegrain the smaller module μ, replacing it with an effective rate ${y}_{\mu}$, and then eliminate all the remaining metabolites in ν. Since effective rates after coarsegraining do not depend on the order of metabolite elimination, ${y}_{\nu}$ must depend on the rate constants ${\overrightarrow{x}}_{\mu}$ only indirectly, through ${y}_{\mu}$. In other words, all the information about the smaller module μ that is relevant for the performance of the larger module ν is contained in ${y}_{\mu}$. 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}_{\mu}$ 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, ${\delta}^{A}{y}_{\mu}$, ${\delta}^{B}{y}_{\mu}$ and $\epsilon {y}_{\mu}$. 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}_{\mu}$ that represents module μ within the larger module ν, and consider ${y}_{\nu}$ as a function of ${y}_{\mu}$, as in Equation 4. To obtain ${f}_{1}({y}_{\mu})$, 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 coarsegraining 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 $f}_{1$ 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 secondorder 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 $\epsilon {y}_{\mu}$ will occur, we need to obtain ${y}_{\mu}$ 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 coarsegraining 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 $f}_{2$ 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 serialparallel (see Figure 3). Second, Proposition 4 and Corollary 4 show that coarsegraining does not alter the strict relationships, that is, if reactions $a$ and $b$ are strictly serial or strictly parallel before coarsegraining they will remain so after coarsegraining. 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 ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$ in Figure 8), then they are always strictly parallel, no matter what the rest of the module looks like. Evaluating the first and secondorder control coefficients for the function $f}_{2$ that corresponds to this topological class reveals that $\epsilon {y}_{\mu}\le 0$ for any parameter values of $f}_{2$.
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 serialparallel (e.g., class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{\varnothing}}$). Consequently, the sign of $\epsilon {y}_{\mu}$ 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 subclasses define parametric subfamilies of function $f}_{2$, and we can explicitly calculate $\epsilon {y}_{\mu}$ for all such functions. However, such bruteforce 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 serialparallel, 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 subfamily of function $f}_{2$, and I explicitly evaluate the first and secondorder control coefficients for all these subfamilies (see Proposition 7 and Proposition 8) which essentially completes the proof of Theorem 2.
Network coarsegraining
Notations and definitions
Request a detailed protocolHere, 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 $i\leftrightarrow j$ has one substrate $i\in A$ and one product $j\in A$, and it is fully described by its rate constant ${x}_{ij}$. By definition, ${x}_{ii}=0$. I denote the set of all reactions by $R=\{i\leftrightarrow j:i,j\in A,{x}_{ij}>0\}$. The dynamics of metabolite concentrations ${S}_{1},\mathrm{\dots},{S}_{n}$ in the network $\mathcal{N}$ are governed by equations
where
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}_{\mu}=\{i\leftrightarrow j\in R:i,j\in {A}_{\mu}\cup {A}_{\mu}^{\mathrm{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}_{\mu}^{\mathrm{i}}\subset {R}_{\mu}$. For example, module μ in Figure 1A has one internal reaction $3\leftrightarrow 4$. 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}_{\mu}^{\mathrm{io}}\subset {R}_{\mu}$. (I reserve uppercase ‘I/O’ for metabolites and use lowercase ‘i/o’ for reactions.) For example, module μ in Figure 1A has three i/o reactions $1\leftrightarrow 3$, $1\leftrightarrow 4$ and $2\leftrightarrow 4$. 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}_{\mu}^{\mathrm{b}}\subset {R}_{\mu}$. For example, module μ in Figure 1A has no bypass reactions but reaction $1\leftrightarrow 5$ is a bypass reaction for module ν. By definition, all these three sets of reactions ${R}_{\mu}^{\mathrm{i}}$, ${R}_{\mu}^{\mathrm{io}}$ and ${R}_{\mu}^{\mathrm{b}}$ are nonoverlapping, and ${R}_{\mu}={R}_{\mu}^{\mathrm{i}}\cup {R}_{\mu}^{\mathrm{io}}\cup {R}_{\mu}^{\mathrm{b}}$.
Another important concept are the simple paths that lie within a module. For any two metabolites $i,j\in A\cup {A}_{\mu}^{\mathrm{IO}}$, I denote a simple path between them that lies within μ as ${p}_{ij}^{\mu}$ or, equivalently as $i\leftrightarrow {k}_{1}\leftrightarrow \mathrm{\dots}\leftrightarrow {k}_{m}\leftrightarrow j$ (where all ${k}_{\mathrm{\ell}}\in {A}_{\mu}$). I say that each of the metabolites ${k}_{\mathrm{\ell}}$ belongs to (or is contained in) path ${p}_{ij}^{\mu}$ (denoted by ${k}_{\mathrm{\ell}}\in {p}_{ij}^{\mu}$). Similarly, I say that each of the reactions ${k}_{\mathrm{\ell}}\leftrightarrow {k}_{\mathrm{\ell}+1}$ belong to (or are contained in) path ${p}_{ij}^{\mu}$ (denoted by ${k}_{\mathrm{\ell}}\leftrightarrow {k}_{\mathrm{\ell}+1}\in {p}_{ij}^{\mu}$). I will drop superindex μ from ${p}_{ij}^{\mu}$ 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 protocolIn this section, I formally introduce and characterize the coarsegraining 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 starmesh transformation in the theory of electric circuits (Versfeld, 1970). The resulting network is a coarsegrained 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 protocolI 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 coarsegrained 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 $\mu =({A}_{\mu},{\overrightarrow{x}}_{\mu})$ contains $m$ internal metabolites. Let $k\in {A}_{\mu}$ be the internal metabolites that will be eliminated. Let ${A}^{\{k\}}=A\setminus \{k\}$ be the reduced metabolite set and let ${\overrightarrow{x}}^{\{k\}}$ be the reduced $(n1)\times (n1)$ matrix of rate constants defined by
where ${D}_{k}$ 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, ${x}_{ij}^{\{k\}}={x}_{ij}$ 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 ${K}_{ij}={K}_{i\mathrm{\ell}}{K}_{\mathrm{\ell}j}$ for any metabolites $i,j,\mathrm{\ell}$, the rate constants ${x}_{ij}^{\{k\}}$ obey Haldane’s relationships. Therefore, the reduced metabolite set ${A}^{\{k\}}$ and the reduced rate matrix ${\overrightarrow{x}}^{\{k\}}$ define a new ‘coarsegrained’ metabolic network ${\mathcal{N}}^{\{k\}}=({A}^{\{k\}},{\overrightarrow{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 coarsegrained network ${\mathcal{N}}^{\{k\}}$ depends only on the reaction set of this module in the original network $\mathcal{N}$, but not on the particular values in the rate matrix ${\overrightarrow{x}}_{\mu}$.
Next, I will show that the dynamics of metabolites in the coarsegrained network ${\mathcal{N}}^{\{k\}}$ are identical to the dynamics of metabolites in the original network $\mathcal{N}$ where metabolite $k$ is at QSS. Note that if metabolite $k$ is at QSS in the network $\mathcal{N}$, its concentration is given by
which follows from Equation 10. Now, the dynamics of metabolites in the network ${\mathcal{N}}^{\{k\}}$ are governed by equations
where ${D}_{i}^{\{k\}}={\sum}_{j\in {A}^{\{k\}}}{x}_{ij}^{\{k\}}$. As mentioned above, ${x}_{ij}^{\{k\}}={x}_{ij}$ 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 $\mathcal{N}$. Next, consider the dynamics of the I/O and internal metabolites for module μ in the coarsegrained network ${\mathcal{N}}^{\{k\}}$, that is, those in the set ${A}_{\mu}^{\mathrm{IO}}\cup {A}_{\mu}\setminus \{k\}$. For any such metabolite $i$, the sum in the righthand side of Equation 15 can be rewritten as
According to Equation 14, the second term in Equation 16 equals ${x}_{ki}{S}_{k}$, so that Equation 16 becomes
For any metabolite $i\in {A}_{\mu}^{\mathrm{IO}}\cup {A}_{\mu}\setminus \{k\}$, the second term in the righthand side of Equation 15 can be rewritten as
Substituting Equation 17 and Equation 18 into Equation 15, we see that Equation 15 is in fact equivalent to Equation 10 for all $i\in A\setminus \{k\}$ with ${S}_{k}$ given by Equation 14.
The coarsegraining procedure (CGP)
Here, I define the CGP for an arbitrary set of internal metabolites by applying metabolite elimination recursively.
Let $E\subseteq {A}_{\mu}$ be an arbitrary subset of metabolites internal to module μ in the metabolic network $\mathcal{N}$ and let ${n}_{E}$ be the number of elements in $E$. Let ${A}^{E}=A\setminus E$ be the reduced metabolite set after the metabolites have been eliminated. I define the reduced $(n{n}_{E})\times (n{n}_{E})$ matrix of rate constants ${\overrightarrow{x}}^{E}$ as follows. If ${n}_{E}=1$, the matrix ${\overrightarrow{x}}^{E}$ is defined by Equation 12 and Equation 13. If ${n}_{E}>1$, then I define it recursively. Suppose that all metabolites in $E$ other than some metabolite $k\in E$ have been previously eliminated, such that the coarsegrained network ${\mathcal{N}}^{{E}^{\prime}}=({A}^{{E}^{\prime}},{\overrightarrow{x}}^{{E}^{\prime}})$ is defined, with the set of eliminated metabolites ${E}^{\prime}=E\setminus \{k\}$, ${A}^{{E}^{\prime}}=A\setminus {E}^{\prime}$, and the known matrix ${\overrightarrow{x}}^{{E}^{\prime}}$. Then, I define the matrix ${\overrightarrow{x}}^{E}$ through the elimination of metabolite $k$ from ${\mathcal{N}}^{{E}^{\prime}}$, that is,
with
I define the the coarsegraining procedure that eliminates the metabolite set E as a map
As with the elimination of a single metabolite, it is straightforward to show that the rate constants ${x}_{ij}^{E}$ obey Haldane’s relationships, so that ${\mathcal{N}}^{E}$ is indeed a metabolic network. ${\mathrm{CG}}^{E}$ maps module μ within the metabolic network $\mathcal{N}$ onto a subnetwork ${\mu}^{\prime}$ within the metabolic network ${\mathcal{N}}^{E}$. It is straightforward to show that ${\mu}^{\prime}$ is a module. Whenever there is no ambiguity, I will label both the original and the coarsegrained versions of the module by μ. To simplify notations, if the CGP eliminates the entire module μ (i.e., if $E={A}_{\mu}$), I label it ${\mathrm{CG}}^{\mu}$. I label the coarsegrained network that restults from the application of ${\mathrm{CG}}^{\mu}$ by ${\mathcal{N}}^{\mu}$ and I label the effective rate of the reaction substituting module μ in network ${\mathcal{N}}^{\mu}$ as ${y}_{\mu}$.
Intuitively, the result of coarsegraining 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 ${x}_{ij}^{E}$. First, by applying the recursion Equation 19, it is easy to show that elimination of two metabolites $E=\{k,\mathrm{\ell}\}$ yields effective rate constants
As expected, Equation 22 and Equation 23 are symmetric with respect to the eliminated metabolites k and $\mathrm{\ell}$. Extrapolating from Equation 23, it is possible to show that for an arbitrary metabolite subset $E\subseteq {A}_{\mu}$ that contains ${n}_{E}$ metabolites,
Here, the second sum is taken over all ${n}_{E}!/({n}_{E}L)!$ ordered lists of metabolites $({k}_{1},\mathrm{\dots},{k}_{L})$ 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 ${x}_{ij}^{E}$ 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 coarsegraining procedure.
The coarsegraining procedure ${\mathrm{CG}}^{E}$ eliminating metabolites in set $E\subseteq {A}_{\mu}$ has the following useful properties which follow from Equation 24.
The effective rate constants ${x}_{ij}^{E}$ do not depend on the order in which metabolites are eliminated. Therefore, the composition of coarsegraining procedures is commutative, that is, if $E={E}_{1}\cup {E}_{2}$, where $E}_{1$ and $E}_{2$ are two nonoverlapping subsets of ${A}_{\mu}$, then
${\mathrm{CG}}^{{E}_{1}}\circ {\mathrm{CG}}^{{E}_{2}}={\mathrm{CG}}^{{E}_{2}}\circ {\mathrm{CG}}^{{E}_{1}}={\mathrm{CG}}^{E}.$If at least one of the metabolites $i$ or $j$ is not adjacent to any of the eliminated metabolites, then ${x}_{ij}^{E}={x}_{ij}$. In particular, ${x}_{ij}^{E}={x}_{ij}$ if either $i$ and/or $j$ are external to μ.
The topology of module μ after the application of ${\mathrm{CG}}^{E}$ depends only on its original topology but not on the values of its rate constants ${\overrightarrow{x}}_{\mu}$.
If both metabolites $i$ and $j$ are adjacent to at least one eliminated metabolite, then ${x}_{ij}^{E}={x}_{ij}+\alpha ,$ where $\alpha \ge 0$ depends only on the rate constants of reactions that involve at least one eliminated metabolite. In particular, if both $k$ and $\mathrm{\ell}$ are from $A\setminus E$, then ${x}_{ij}^{E}$ is independent of ${x}_{k\mathrm{\ell}}$.
If $E={A}_{\mu}$, then the effective rate constant ${y}_{\mu}$ of module μ depends on the rate matrix ${\overrightarrow{x}}_{\mu}$ but does not depend on any other reaction rates, that is, Equation 1 holds. Furthermore, the functional form of $F({\overrightarrow{x}}_{\mu})$ depends only on the topology of module μ, that is, all modules with the same topology are mapped onto ${y}_{\mu}$ by the same function $F$.
Suppose that module μ is nested in a larger module $\nu =({A}_{\nu},{\overrightarrow{x}}_{\nu})$ (see Figure 1A). It follows from Property #1 that ${\mathrm{CG}}^{\nu}={\mathrm{CG}}^{\mu}\circ {\mathrm{CG}}^{{A}_{\nu}\setminus {A}_{\mu}}$, that is, ${y}_{\nu}$ 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}_{\mu}$ and then eliminating the remaining metabolites in ${A}_{\nu}$. In the network ${\mathcal{N}}^{\mu}$ after the first stage of coarsegraining, all rate constants ${\overrightarrow{x}}_{\nu \setminus \mu}$ are identical to those in the original network, that is, they are independent of ${\overrightarrow{x}}_{\mu}$, by virtue of Property #2. Therefore ${y}_{\nu}$ depends on the rate constants ${\overrightarrow{x}}_{\mu}$ of reactions within module μ only through the effective rate constant ${y}_{\mu}$ of module μ. In other words, Equation 4 holds.
If, in the original network $\mathcal{N}$, 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 coarsegrained network ${\mathcal{N}}^{E}$.
If, in the metabolic network $\mathcal{N}$, metabolites $i$ and $j$ are not adjacent (i.e. ${x}_{ij}=0$) and no simple path exists within the set $E$ (i.e. such that all nonterminal metabolites in this path are from $E$) that connects metabolites $i$ and $j$, then metabolites $i$ and $j$ are also not adjacent in the coarsegrained network ${\mathcal{N}}^{E}$ (i.e., ${x}_{ij}^{E}=0$).
It follows from properties #7 and #8 that for a simple path ${p}_{{\mathrm{\ell}}_{1}{\mathrm{\ell}}_{m}}={\mathrm{\ell}}_{1}\leftrightarrow {\mathrm{\ell}}_{2}\leftrightarrow \mathrm{\cdots}\leftrightarrow {\mathrm{\ell}}_{m}$ to exist in module μ after the application of ${\mathrm{CG}}^{\mu}$, it is neccessary and sufficient that for each $i=1,\mathrm{\dots},m1$, either ${\mathrm{\ell}}_{i}$ and ${\mathrm{\ell}}_{i+1}$ are adjacent in the original module μ or there exists a simple path ${\mathrm{\ell}}_{i}\leftrightarrow \mathrm{\cdots}\leftrightarrow {\mathrm{\ell}}_{i+1}$ within the original module μ all of whose nonterminal 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 coarsegraining procedures, as follows. Suppose that ${\mathrm{CG}}^{{E}_{1}}$ and ${\mathrm{CG}}^{{E}_{2}}$ are two coarsegraining procedures of network $\mathcal{N}$ for two subsets of metabolites ${E}_{1}\subset {A}_{\mu}$ and ${E}_{2}\subset {A}_{\mu}$. If the sets $E}_{1$ and $E}_{2$ are nonoverlapping, ${\mathrm{CG}}^{{E}_{2}}$ is also defined for the coarsegrained network ${\mathcal{N}}^{{E}_{1}}$ which is the result of applying ${\mathrm{CG}}^{{E}_{1}}$ to the original network $\mathcal{N}$. The result of applying ${\mathrm{CG}}^{{E}_{2}}$ to the ${\mathcal{N}}^{{E}_{1}}$ is called the composition of coarsegraining procedures ${\mathrm{CG}}^{{E}_{\mathrm{1}}}$ and ${\mathrm{CG}}^{{E}_{2}}$ of the original network $\mathcal{N}$ and is denoted as ${\mathrm{CG}}^{{E}_{1}}\circ {\mathrm{CG}}^{{E}_{2}}$.
As defined above, coarsegraining 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 steadystate 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 protocolLet $E$ be any subset of metabolites internal to module μ. Then,
There exists a joint QSS solution ${\overline{S}}_{i}$ for all metabolites $i\in E$, given the concentrations of the remaining internal and I/O metabolites for module μ.
The dynamics of all remaining metabolites in $A\setminus E$ in the coarsegrained metabolic network ${\mathcal{N}}^{E}$ are the same as in $\mathcal{N}$ where all metabolites in $E$ are at QSS.
Corollary 1
Request a detailed protocolWithout loss of generality, suppose that the I/O metabolites for module μ are labeled 1 and 2 and its internal metabolites are labeled ${A}_{\mu}\mathrm{=}\mathrm{\{}\mathrm{3}\mathrm{,}\mathrm{4}\mathrm{,}\mathrm{\dots}\mathrm{,}m\mathrm{\}}$. There exists a unique QSS ${\overline{S}}_{i}$ for all $i\mathrm{\in}{A}_{\mu}$. The QSS concentrations can be obtained by recursively applying equation.
for $k\mathrm{=}\mathrm{3}\mathrm{,}\mathrm{4}\mathrm{,}\mathrm{\dots}\mathrm{,}m$.
Equation 25 follows from Equation 10 for the coarsegrained network obtained by eliminating metabolites $k+1,\mathrm{\dots},m$.
Corollary 2
Request a detailed protocolWithout 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}_{\mu}$ can be calculated using Equation 19 and Equation 20 or Equation 24. The dynamics of all metabolites in the resulting coarsegrained metabolic network are identical to their dynamics in the original network $\mathrm{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 protocolCorollary 2 provides a method for replacing any module μ at QSS with an effective rate ${y}_{\mu}=F({\overrightarrow{x}}_{\mu})$, 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 protocolConsider a linear pathway with I/O metabolites 1 and $m$ and internal metabolites $2,3,\mathrm{\dots},m1$ (Figure 6A). This labeling of metabolites is more convenient for the linear pathway. To calculate ${y}_{\mu}$, I will apply recursion Equation 19 and Equation 20. I start by eliminating metabolite 2. After this initial coarsegraining step, the resulting module is still a linear pathway, where two reactions $1\leftrightarrow 2\leftrightarrow 3$ were replaced with a single reaction $1\leftrightarrow 3$ with the effective rate constant.
All other rate constants remain unchanged. Next, I eliminate metabolite 3. The resulting module is still a linear pathway, where now three reactions $1\leftrightarrow 2\leftrightarrow 3\leftrightarrow 4$ were replaced with a single reaction $1\leftrightarrow 4$ with the effective rate constant
All other rate constants remain unchanged. Continuing this process until all internal metabolites are eliminated, I obtain
which is identical to the expression originally obtained by Kacser and Burns, 1973.
Two parallel pathways
Request a detailed protocolConsider 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$, $\mathrm{\ell}=4$. Since ${x}_{12}={x}_{34}=0$, Equation 22 simplifies to
Thus, the contributions of parallel pathways are simply added.
Module μ in Figure 1
Request a detailed protocolTo obtain the effective rate constant for module μ shown in Figure 1, I again use Equation 22 with $i=1$, $j=2$, $k=3$, $\mathrm{\ell}=4$.
with ${D}_{3}={x}_{31}+{x}_{32}+{x}_{34}$ and ${D}_{4}={x}_{41}+{x}_{42}+{x}_{43}$.
Classification of modules with respect to ‘marked’ reactions, and the parametric families of functions f_{1} and f_{2}
The CGP described above allows us to calculate the function $F$ that maps the rate matrix ${\overrightarrow{x}}_{\mu}$ for an arbitrary module μ onto the module’s effective rate constant ${y}_{\mu}$. $F$ is a multivariate function of the entire matrix ${\overrightarrow{x}}_{\mu}$. 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}_{\mu}$ depends on these one or two perturbed reactions. I refer to such singledout reactions as ‘marked’. When ${y}_{\mu}$ is considered as a function of the rate constant ξ of one marked reaction, I write ${y}_{\mu}={f}_{1}(\xi )$, as in Equation 2. When ${y}_{\mu}$ is considered as a function of the rate constants ξ and η of two marked reactions, I write ${y}_{\mu}={f}_{2}(\xi ,\eta )$ as in Equation 3.
The functional form of $F$ and, as a consequence, the functional forms of $f}_{1$ and $f}_{2$ 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 $f}_{1$ and $f}_{2$, such that each topology of module μ defines a parametric family of functions $f}_{1$ and $f}_{2$, 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 $f}_{1$ and $f}_{2$ 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 coarsegrained 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 $f}_{1$. The topologies of the three minimal modules are sufficiently simple that the three corresponding parametric functional forms of $f}_{1$ can be easily computed by applying the coarsegraining Equation 19 or Equation 22. This result is formulated in Proposition 2.
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 $f}_{2$ from the same parametric family. These families are characterized in Corollary 3.
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 $\mu =({A}_{\mu},{\overrightarrow{x}}_{\mu})$ and let $a={i}_{a}\leftrightarrow {j}_{a}$ and $b={i}_{b}\leftrightarrow {j}_{b}$ be two reactions from its set of reactions ${R}_{\mu}$. I call a pair $(\mu ,a)$ a singlemarked module and I call a triplet $(\mu ,a,b)$ a doublemarked module, and I refer to reactions $a$ and $b$ as marked within module μ. The topology of a singlemarked module $(\mu ,a)$ is determined not only by the reaction matrix ${R}_{\mu}$, but also by the position of the marked reaction, so I refer to the pair $({R}_{\mu},a)$ as the topology of the singlemarked module $\mathrm{(}\mu \mathrm{,}a\mathrm{)}$. Similarly, I refer to the triplet $({R}_{\mu},a,b)$ as the topology of the doublemarked module $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$. I denote by ${\overrightarrow{x}}_{\mu \setminus a}$ the matrix of rate constants of all reactions in module μ other than reaction $a$ and I denote by ${\overrightarrow{x}}_{\mu \setminus \{a,b\}}$ the matrix of all rate constants in module μ other than reactions $a$ and $b$. I denote the sets of all single and doublemarked modules by ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$, respectively. To avoid metabolite labeling ambiguities, I adopt the following conventions:
The I/O metabolites are labeled 1 and 2 and the internal metabolites are labeled $3,4,\mathrm{\dots}$;
${i}_{a},{j}_{a}\in \{1,2,3,4\}$; ${i}_{b},{j}_{b}\in \{1,2,3,4,5,6\}$;
${i}_{a}<{j}_{a}$, ${i}_{b}<{j}_{b}$, ${i}_{a}\le {i}_{b}$, ${j}_{a}\le {j}_{b}$.
It is easy to see that the set of all singlemarked modules ${\mathcal{M}}_{1}$ can be partitioned into three nonoverlapping topological classes depending on the type of the marked reaction $a$. I denote the classes of all singlemarked modules where the marked reaction is bypass, i/o or internal (see Notations and definitions) by ${\mathcal{M}}^{\mathrm{b}}$, ${\mathcal{M}}^{\mathrm{io}}$ and ${\mathcal{M}}^{\mathrm{i}}$, respectively (Figure 7). Similarly, the set ${\mathcal{M}}_{2}$ can be partitioned into nine nonoverlapping 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.
Each topological class contains infinitely many modules, with various numbers of metabolites and various topologies. However, for each topological class $\mathcal{M}$, there is a minimum number of internal metabolites ${m}_{\mathcal{M}}$, such that all modules within $\mathcal{M}$ must have at least ${m}_{\mathcal{M}}$ internal metabolites. I denote the set of metabolites minimal in the topological class $\mathcal{M}$ by ${A}_{\mathcal{M}}$. It is clear that for the singlemarked module classes ${\mathcal{M}}^{\mathrm{b}}$ and for ${\mathcal{M}}^{\mathrm{io}}$, ${m}_{{\mathcal{M}}^{\mathrm{b}}}={m}_{{\mathcal{M}}^{\mathrm{io}}}=1$ and ${A}_{{\mathcal{M}}^{\mathrm{b}}}={A}_{{\mathcal{M}}^{\mathrm{io}}}=\{3\}$, and for ${\mathcal{M}}^{\mathrm{i}}$, ${m}_{{\mathcal{M}}^{\mathrm{i}}}=2$ and ${A}_{{\mathcal{M}}^{\mathrm{i}}}=\{3,4\}$ (see second column in Figure 7). For the doublemarked modules, ${m}_{\mathcal{M}}$ and ${A}_{\mathcal{M}}$ are given in Table 1 and illustrated in Figure 8.
If a singlemarked module from the topological class $\mathcal{M}$ has the minimum number of metabolites ${n}_{\mathcal{M}}$ in that class, I call such module and its topology minimal in $\mathcal{M}$. There may be several minimal topologies in a topological class, but there is only one minimal topology that is fully connected. A topology $({R}_{\mu},a)$ is called fully connected if the reaction set ${R}_{\mu}$ is complete in the sense that it contains reactions between all pairs of metabolites in the minimal metabolite set ${A}_{\mathcal{M}}$. I denote such complete reaction set for the class $\mathcal{M}$ by ${R}_{\mathcal{M}}$, and I denote the respective fully connected topology by $({R}_{\mathcal{M}},a)$. I employ the same terminology and analogous notations for doublemarked 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 $f}_{1$ and characterizes them. The idea of the proof is the following. According to Property 5 (Box 1), all singlemarked modules that are mapped by the CGP onto a minimal module with the same topology $({R}_{\mu},a)$ have the same functional form of $f}_{1$. In other words, each minimal topology $({R}_{\mu},a)$ specifies a parameteric family of the function $f}_{1$. Since the number of possible minimal topologies is finite, the claim of Theorem 1 can be tested for each corresponding functional form of $f}_{1$. However, the number of minimal topologies is rather large. Fortunately, another simplification is possible. Since the reaction set ${R}_{\mu}$ of any minimal singlemarked module is a subset of the complete reaction set ${R}_{\mathcal{M}}$, the fully connected topology $({R}_{\mathcal{M}},a)$ specifies the largest parametric family of the function $f}_{1$ for the class $\mathcal{M}$, 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 singlemarked modules that belong to the topological class $\mathcal{M}$ are described by functions $f}_{1$ that belong to one parameteric family corresponding to the fully connected topology minimal in $\mathcal{M}$. The three parameteric families of $f}_{1$ 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 $f}_{1$ 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 protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{)}$ be a singlemarked module, and let ξ be the rate constant of reaction $a$. Then ${y}_{\mu}\mathrm{=}{f}_{\mathrm{1}}\mathit{}\mathrm{(}u\mathrm{)}$, where $u\mathrm{=}\xi \mathrm{+}\alpha $ for some $\alpha \mathrm{\ge}\mathrm{0}$, and the function ${f}_{\mathrm{1}}$ is given by one of the following expressions.
Here ${D}_{\mathrm{3}}\mathrm{=}{w}_{\mathrm{31}}\mathrm{+}{w}_{\mathrm{32}}\mathrm{+}u$, ${D}_{\mathrm{4}}\mathrm{=}{w}_{\mathrm{41}}\mathrm{+}{w}_{\mathrm{42}}\mathrm{+}u\mathrm{/}{K}_{\mathrm{34}}$, and all ${w}_{i\mathit{}j}\mathrm{\ge}\mathrm{0}$ are independent of ξ.
Proof
Request a detailed protocolSince any singlemarked module $(\mu ,a)$ belongs to exactly one of three classes ${\mathcal{M}}^{\mathrm{b}}$, ${\mathcal{M}}^{\mathrm{io}}$ and ${\mathcal{M}}^{\mathrm{io}}$, I consider these three cases one by one.
Case $\mathrm{(}\mu \mathrm{,}a\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{b}}$. According to the labeling conventions outlined above, $a=1\leftrightarrow 2$ (see Figure 7). Equation 29 follows directly from Property #4 of the CGP (Box 1).
Case $\mathrm{(}\mu \mathrm{,}a\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}}$. According to the labeling conventions, $a=1\leftrightarrow 3$ (see Figure 7). According to Property #1 of the CGP, module μ can be coarsegrained in two stages, by first applying ${\mathrm{CG}}^{{A}_{\mu}\setminus \{3\}}$ which eliminates metabolites $4,\mathrm{\dots},m$ (those that do not participate in the marked reaction) and then applying ${\mathrm{CG}}^{\{3\}}$ which eliminates the remaining metabolite 3. Mathematically, ${\mathrm{CG}}^{\mu}={\mathrm{CG}}^{{A}_{\mu}\setminus \{3\}}\circ {\mathrm{CG}}^{\{3\}}$. After applying ${\mathrm{CG}}^{{A}_{\mu}\setminus \{3\}}$, the resulting coarsegrained module ${\mu}^{\prime}$ has a single internal metabolite 3 and at most three effective reactions $1\leftrightarrow 2$, $1\leftrightarrow 3$ and $2\leftrightarrow 3$ (Figure 7), that is, it is minimal in ${\mathcal{M}}^{\mathrm{io}}$. By virtue of Properties #2 and #4 of the CGP, the effective rate constants $w}_{12$, $w}_{23$ are independent of ξ and $u\equiv {w}_{13}=\xi +\alpha $. Note that $w}_{12$ may equal zero, but ${w}_{23}\ne 0$ because ${\mu}^{\prime}$ is a module. Regardless, the reaction set ${R}_{{\mu}^{\prime}}$ of module ${\mu}^{\prime}$ is always a subset of the complete reaction set ${R}_{{\mathcal{M}}^{\mathrm{io}}}$. Thus, to obtain the effective rate constant ${y}_{\mu}$, I consider the most general case when ${\mu}^{\prime}$ is fully connected and eliminate the remaining internal metabolite 3, which leads to Equation 30.
Case $(\mu ,a)\in {\mathcal{\mathcal{M}}}^{\mathrm{i}}$. According to the labeling conventions, $a=3\leftrightarrow 4$ (see Figure 7). Otherwise, the logic of the proof is exactly the same as for the case $(\mu ,a)\in {\mathcal{M}}^{\mathrm{io}}$.
Next I prove Proposition 3 which states that, for any doublemarked module that belongs to a given topological class, there exists a doublemarked module that is minimal in the same class, such that both modules have the same function $f}_{2$. The corresponding minimal module is obtained from the original module by applying the CGP. This proposition is important because it implies that all functions $f}_{2$ can be completely characterized by only examing minimal modules. Then, analogously to singlemarked modules, Corollary 3 states that function $f}_{2$ 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 protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ be a doublemarked module that belongs to the topological class $\mathrm{M}$, and let ξ and η be the rate constants of reactions $a$ and $b$, respectively. Then there exist nonnegative constants α and β and a doublemarked module $\mathrm{(}{\mu}^{\mathrm{\prime}}\mathrm{,}a\mathrm{,}b\mathrm{)}$ minimal in $\mathrm{M}$ such that ${y}_{\mu}\mathrm{=}{y}_{{\mu}^{\mathrm{\prime}}}\mathrm{=}{f}_{\mathrm{2}}\mathit{}\mathrm{(}u\mathrm{,}v\mathrm{)}$, where
are the rate constants of the marked reactions $a$ and $b$ in ${\mu}^{\mathrm{\prime}}$, respectively, and all other rate constants in ${\mu}^{\mathrm{\prime}}$ are independent of ξ or η. Module ${\mu}^{\mathrm{\prime}}$ is obtained from μ by the coarsegraining procedure ${\mathrm{CG}}^{\mu \mathrm{\setminus}\mathrm{\{}a\mathrm{,}b\mathrm{\}}}$ that eliminates all metabolites internal to module μ that do not participate in reactions $a$ or $b$.
Proof
Request a detailed protocolTo prove this proposition, I will construct the doublemarked module $({\mu}^{\prime},a,b)$ minimal in $\mathcal{M}$ by applying ${\mathrm{CG}}^{\mu \setminus \{a,b\}}$. Let ${m}_{\mathcal{M}}$ be the mimimal number of internal metabolites in class $\mathcal{M}$ (see Table 1). According to the metabolite labeling conventions, metabolites ${n}_{\mathcal{M}}+3,{n}_{\mathcal{M}}+4,\mathrm{\dots}$ are neither I/O nor do they participate in the marked reactions. ${\mathrm{CG}}^{\mu \setminus \{a,b\}}$ eliminates all these metabolites and maps module μ onto module ${\mu}^{\prime}$, all of whose internal metabolites participate in reactions $a$ and/or $b$. Therefore, $({\mu}^{\prime},a,b)$ is minimal in class $\mathcal{M}$ (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 ${\mu}^{\prime}$ 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}_{\mu}={y}_{{\mu}^{\prime}}$ follows from Property #1 of the CGP, ${\mathrm{CG}}^{\mu}={\mathrm{CG}}^{\mu \setminus \{a,b\}}\circ {\mathrm{CG}}^{{A}_{{\mu}^{\prime}}}$.
Corollary 3
Request a detailed protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ be a doublemarked module, and let ξ and η be the rate constants of reactions $a$ and $b$, respectively. The function ${f}_{\mathrm{2}}$ that maps ξ and η onto module’s effective rate constant ${y}_{\mu}$ belongs to one of nine parametric families. If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{b}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{IO}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{b}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{I}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{IO}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{\varnothing}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}}$, then
If $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}}$, then
In Equation 34 through Equation 35, $u$ and $v$ are given by Equation 32 and Equation 33. All effective activities ${w}_{i\mathit{}j}\mathrm{\ge}\mathrm{0}$ are constants (other than cases where they stand for $u$ or $v$) that depend on the topology of module μ and on the parameters ${\overrightarrow{x}}_{\mu \mathrm{\setminus}\mathrm{\{}a\mathrm{,}b\mathrm{\}}}$ but do not depend on ξ and η.
Proof
Request a detailed protocolThis statement and Equation 34 through Equation 35 follow directly from Proposition 3 and the fact that the reaction set of any doublemarked 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 protocolConsider a higherlevel phenotype $y$, such as the effective activity of a module, which is function of a multivariate lowerlevel phenotype $\overrightarrow{x}=({x}_{1},{x}_{2},\mathrm{\dots},{x}_{n})$, such as the rates of individual reactions within the module, $y=F\left(\overrightarrow{x}\right)$. Denote the wildtype values of the phenotypes as ${\overrightarrow{x}}^{0}=({x}_{1}^{0},{x}_{2}^{0},\mathrm{\dots},{x}_{n}^{0})$ and ${y}^{0}=F\left({\overrightarrow{x}}^{0}\right)$. Consider a mutation that perturbes these values, so that the mutant has lowerlevel phenotypic values ${\overrightarrow{x}}^{\prime}=({x}_{1}^{\prime},{x}_{2}^{\prime},\mathrm{\dots},{x}_{n}^{\prime})$. The relative effect of the mutation on phenotype x_{i} is $\delta {x}_{i}={x}_{i}^{\prime}/{x}_{i}^{0}1$. If all $\parallel \delta {x}_{i}\parallel \ll 1$ where $\parallel \overrightarrow{x}\parallel $ denotes the length of vector $\overrightarrow{x}$, then the value of the higherlevel phenotype ${y}^{\prime}$ in the mutant is given by
where
which I refer to as first and secondorder control coefficients of the lowerlevel phenotypes $x}_{i$ and $x}_{j$ with respect to the higherlevel phenotype $y$.
Now consider two single mutants, A and B, and the doublemutant AB. Each mutation A and B and their combination may perturb all $x}_{i$ phenotypes such that ${x}_{i}^{A}={x}_{i}^{0}\left(1+{\delta}^{A}{x}_{i}\right)$, ${x}_{i}^{B}={x}_{i}^{0}\left(1+{\delta}^{B}{x}_{i}\right)$, and ${x}_{i}^{AB}={x}_{i}^{0}\left(1+{\delta}^{AB}{x}_{i}\right)={x}_{i}^{0}\left(1+{\delta}^{A}{x}_{i}+{\delta}^{B}{x}_{i}+2{\delta}^{A}{x}_{i}{\delta}^{B}{x}_{i}\epsilon {x}_{i}\right).$
Assuming that $\parallel {\delta}^{A}\overrightarrow{x}\parallel \ll 1$, $\parallel {\delta}^{B}\overrightarrow{x}\parallel \ll 1$ and $\parallel {\delta}^{AB}\overrightarrow{x}\parallel \ll 1$, using the approximation in Equation 44 and the definition of $\epsilon {x}_{i}$ (analogous to Equation 5), I obtain
where $o(1)$ refers to terms that are vanishingly small as ${\delta}^{A}{x}_{i}\to 0$ , ${\delta}^{B}{x}_{i}\to 0$, $i=1,\mathrm{\dots}n$.
I examine two special cases of Equation 49. The first special case is when both mutations affect a single phenotype $x}_{k$, that is, when all ${\delta}^{A}{x}_{i}=0$ and all ${\delta}^{B}{x}_{i}=0$ except for $i=k$. Then Equation 47, Equation 48, Equation 49 simplify to
Equation 52 is equivalent to Equation 6.
The second special case is when mutation A affects a single phenotypes $x}_{k$ and mutation B affects a single phenotype ${x}_{\mathrm{\ell}}$ ($k\ne \mathrm{\ell}$), i.e., all ${\delta}^{A}{x}_{i}=0$ except for $i=k$, all ${\delta}^{B}{x}_{i}=0$ except for $i=\mathrm{\ell}$, and all $\epsilon {x}_{i}=0$. Then Equation 47, Equation 48, Equation 49 simplify to
Equation 55 is equivalent to Equation 9.
Calculation of epistasis in simple modules
Request a detailed protocolEquation 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 protocolFirst, consider how epistasis propagates through a linear pathway (Figure 6A). For simplicity, assume that both mutations A and B affect the same reaction $1\leftrightarrow 2$. It follows from Equation 26 that
where α is a positive constant. Therefore, the first and secondorder control coefficients of reaction $1\leftrightarrow 2$ with respect to the flux through the linear pathway μ are given by
Substituting these expressions into the expression for the fixed point $\overline{\epsilon}=H{\left(2C\left(1C\right)\right)}^{1}$, I find that $\overline{\epsilon}=1$, irrespectively of the rates of other reactions in the linear pathway. This implies that epistasis $\epsilon \xi <1$ at the level of reaction $1\leftrightarrow 2$ would induce a lower value of epistasis $\epsilon {y}_{\mu}<\epsilon \xi <1$ at the level of the entire linear pathway, any value $\epsilon \xi >1$ would induce a higher value of epistasis $\epsilon {y}_{\mu}>\epsilon \xi >1$, and $\epsilon \xi =1$ would induce $\epsilon {y}_{\mu}=1$.
Now consider emergence of epistasis in a linear pathway. Suppose that mutation A affects reaction $1\leftrightarrow 2$ and mutation B affects reaction $2\leftrightarrow 3$. Denote the rate constant of reactions $1\leftrightarrow 2$ and $2\leftrightarrow 3$ by $\xi \equiv {x}_{12}$ and $\eta \equiv {x}_{23}$, respectively. It follows from Equation 26 that
where β is a positive constant. Therefore,
which, after substituting into Equation 9, yield $\epsilon {y}_{\mu}=1$. Together with the fact that $\overline{\epsilon}=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 protocolSuppose that mutation A affects reaction $1\leftrightarrow 3$ and mutation B affects reaction $1\leftrightarrow 4$ in the linear metabolic pathway shown in Figure 6B. Denote the rate constants of reaction $1\leftrightarrow 3$ and $1\leftrightarrow 4$ by $\xi \equiv {x}_{13}$ and $\eta \equiv {x}_{14}$. It follows from Equation 27 that
where $\alpha =1/({K}_{13}{x}_{32})$ and $\beta =1/({K}_{14}{x}_{42})$. Thus, we have ${H}_{\xi \eta}=0$, and there is no epistasis between such mutations.
Module ν in Figure 3A
Request a detailed protocolI denote the rate of the reactions affected by mutations A and B by $\xi ={x}_{13}$ and $\eta ={x}_{42}$, and I also denote $z={x}_{34}$. I will calculate the epistasis coefficient $\epsilon {y}_{\nu}$ in two stages, by first calculating the epistasis coefficient $\epsilon {y}_{\mu}$ and then propagating it to $\epsilon {y}_{\nu}$ using Equation 6. Here I am specifically interested in how $\epsilon {y}_{\nu}$ depends on the rate constant $z$.
To compute epistasis between mutations A and B at the level of module μ, I rewrite Equation 28 as
where $a={x}_{14}/{K}_{13}+{x}_{32}+z$, ${b}_{\xi}={x}_{32}\left({x}_{41}+z/{K}_{34}\right)$, ${b}_{\eta}={x}_{14}\left({x}_{32}+z\right)$, $c={x}_{14}{x}_{32}z/{K}_{34}$, $d=1/{K}_{13}$, ${e}_{\xi}=\left({x}_{41}+z/{K}_{34}\right)/{K}_{13}$, ${e}_{\eta}={x}_{32}+z$, $f={x}_{32}z/{K}_{34}+{x}_{41}z+{x}_{32}{x}_{41}$. I obtain the following expressions for the first and secondorder control coefficients.
where $D=d\xi \eta +{e}_{\xi}\xi +{e}_{\eta}\eta +f$, ${\stackrel{~}{c}}_{1}={x}_{23}/{K}_{24}+\eta $, ${\stackrel{~}{d}}_{1}={x}_{32}\left({x}_{41}+\eta \right)$, ${\stackrel{~}{c}}_{2}=\xi +{x}_{14}$, ${\stackrel{~}{d}}_{2}={x}_{14}\left(\xi /{K}_{13}+{x}_{32}\right)$. Substituting Equation 56 through Equation 58 into Equation 53 through Equation 55, I obtain
where $\stackrel{~}{a}={\stackrel{~}{c}}_{1}{\stackrel{~}{c}}_{2}$ and $\stackrel{~}{b}=\left(\xi /{K}_{13}+{x}_{32}\right){x}_{14}\eta +\left({x}_{41}+\eta \right){x}_{32}\xi $.
To obtain the expression for $\epsilon {y}_{\nu}$, I coarsegrain module ν by eliminating the only remaining internal metabolite 2 and obtain
I then apply equation Equation 6 with
Figure 3B illustrates how $\epsilon {y}_{\nu}$ changes as a function of $z$. It was generated using the following matrix of rate constants:
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.
$z=0$. According to Equation 61, $\epsilon {y}_{\mu}=0$ and hence $\epsilon {y}_{\nu}=H/2{C}^{2}\le 0$, with $C$ and $H$ given by Equation 62 and Equation 63. It is easy to see that in this case the reactions $1\leftrightarrow 3$ and $2\leftrightarrow 4$ that are affected by mutations are strictly parallel because there is a simple path $1\leftrightarrow 3\leftrightarrow 2\leftrightarrow 5$ that contains only reaction $1\leftrightarrow 3$ and there is a simple path $1\leftrightarrow 4\leftrightarrow 2\leftrightarrow 5$ that contains only reaction $2\leftrightarrow 4$ (Figure 3C).
${x}_{15}={x}_{32}={x}_{14}=0$. In this case, module ν is a linear pathway. Therefore, $\epsilon {y}_{\nu}=1$ independently of $z$, as shown above. This fact can also be obtained directly from Equation 61, Equation 62, Equation 63 and Equation 6.
$z\to \mathrm{\infty}$. 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 $3\leftrightarrow 4$. Thus, it follows from Equation 61 that $\epsilon {y}_{\mu}\to \stackrel{~}{a}/\left({\stackrel{~}{c}}_{1}{\stackrel{~}{c}}_{2}\right)=1$, as expected. Then, according to Theorem 1, $\epsilon {y}_{\nu}\ge 1$.
Proof of Theorem 1
As discussed above, the key step toward the proof is Proposition 2, which states that the function $f}_{1$ 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 protocolLet $a$ be the effective reaction within higherlevel module ν that represents the lowerlevel module μ. To simplify notations, I denote ${y}_{\mu}\equiv \xi $. According to Proposition 2, the functional from of $f}_{1$ depends only on the topological class of the singlemarked module $(\nu ,a)$. So, I consider the three classes one by one.
Case $\mathrm{(}\nu \mathrm{,}a\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{b}}$. From Equation 29, $C=\xi /{y}_{\nu}$ and $H=0$. Therefore, inequalities in Equation 7 and Equation 8 hold.
Case $\mathrm{(}\nu \mathrm{,}a\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{io}}$. From Equation 30,
where $D=(\xi +\alpha )/{K}_{13}+{w}_{32}$. From Equation 64, it is clear that $C\ge 0$. Rewriting Equation 64 as
it is also clear that $C\le 1$ since both ratios in this expression do not exceed 1. From Equation 65 and the fact that $0\le C\le 1$, it follows that $\overline{\epsilon}\ge 0$. To show that $\overline{\epsilon}\le 1$, note that
Therefore,
Therefore, inequalities in Equation 7 and Equation 8 hold.
Case $\mathrm{(}\nu \mathrm{,}a\mathrm{)}\mathrm{\in}{\mathrm{M}}^{\mathrm{i}}$. I rewrite Equation 31 as
with $u=\xi +\alpha $, $D=\stackrel{~}{C}u+\stackrel{~}{D}$, $\stackrel{~}{A}=\left({w}_{13}+{w}_{14}\right)\left({w}_{42}+{w}_{32}/{K}_{34}\right)$, $\stackrel{~}{B}=({w}_{31}+{w}_{32}){w}_{14}{w}_{42}+({w}_{41}+{w}_{42}){w}_{13}{w}_{32}$, $\stackrel{~}{C}=({w}_{31}+{w}_{32})/{K}_{34}+({w}_{41}+{w}_{42})$, $\stackrel{~}{D}=({w}_{31}+{w}_{32})({w}_{41}+{w}_{42})$, which yields
Next, it is straightforward to show that $\stackrel{~}{A}\stackrel{~}{D}\stackrel{~}{B}\stackrel{~}{C}={\left({w}_{41}{w}_{32}{w}_{31}{w}_{42}\right)}^{2}/{K}_{31}\ge 0$, which implies that $C\ge 0$. To show that $C\le 1$, I expand the denominator in Equation 66 and obtain
Therefore, numerator in Equation 66 cannot exceed the denominator. The fact that $\overline{\epsilon}\ge 0$ follows directly from Equation 67 together with $C\le 1$. To show that $\overline{\epsilon}\le 1$, first note that
Therefore,
Hence,
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 serialparallel. Then, Proposition 4 and its Corollary 4 establish that strictly parallel (serial) reactions in $(\mu ,a,b)$ are also strictly parallel (serial) in a minimal module $({\mu}^{\prime},a,b)$, which is obtained from μ by eliminating all metabolites that do not participate in the marked reactions. Next, recall that in both modules $(\mu ,a,b)$ and $({\mu}^{\prime},a,b)$ the same function $f}_{2$ 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 $f}_{2$ 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 $f}_{2$.
Unfortunately, the number of minimal topologies is very large, so that such bruteforce 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 $\epsilon {y}_{\mu}\le 0$ for any minimal module μ with any strictly parallel generating topology. Second, Proposition 8 shows that that $\epsilon {y}_{\mu}\ge 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 protocolConsider 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 ${p}_{ij}^{\mu}=i\leftrightarrow k\leftrightarrow \mathrm{\dots}\leftrightarrow \mathrm{\ell}\leftrightarrow j$. If such path contains reactions ${a}_{1},{a}_{2},\mathrm{\dots}$ and does not contain reactions ${b}_{1},{b}_{2},\mathrm{\dots}$, I denote it as ${p}_{ij}^{\mu}({a}_{1},{a}_{2},\mathrm{\dots},{\overline{b}}_{1},{\overline{b}}_{2},\mathrm{\dots})$. I denote the set of all paths ${p}_{ij}^{\mu}({a}_{1},{a}_{2},\mathrm{\dots},{\overline{b}}_{1},{\overline{b}}_{2},\mathrm{\dots})$ by ${\mathcal{P}}_{ij}^{\mu}({a}_{1},{a}_{2},\mathrm{\dots},{\overline{b}}_{1},{\overline{b}}_{2},\mathrm{\dots})$.
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, ${\mathcal{P}}_{12}^{\mu}(a)\ne \mathrm{\varnothing}$ for any reaction $a\in {R}_{\mu}$. 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 doublemarked module $(\mu ,a,b)$ based on the toplogocial relationship between its marked reactions. This classification is orthogonal to that given in Table 1.
Two reactions $a\in {R}_{\mu}$ and $b\in {R}_{\mu}$ are called parallel within module μ if there exists a simple path ${p}_{12}^{\mu}(a,\overline{b})$ and a simple path ${p}_{12}^{\mu}(b,\overline{a})$. Two reactions $a\in {R}_{\mu}$ and $b\in {R}_{\mu}$ are called serial within module μ if there exist at least one simple path ${p}_{12}^{\mu}(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 serialparallel 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 serialparallel. 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 serialparallel within module μ, I also refer to the doublemarked module $(\mu ,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}_{\mu},a,b)$ as serial, parallel, etc.
Recall that coarsegraining procedure ${\mathrm{CG}}^{\mu \setminus \{a,b\}}$ eliminates all metabolites internal to module μ other than those participating in reactions $a$ and $b$. If the doublemarked module $(\mu ,a,b)$ belongs to the topological class $\mathcal{M}$, then, according to Proposition 3, ${\mathrm{CG}}^{\mu \setminus \{a,b\}}$ maps $(\mu ,a,b)$ onto a minimal doublemarked module $({\mu}^{\prime},a,b)$ from the same class $\mathcal{M}$. 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 protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ be a doublemarked module from the topological class $\mathrm{M}$ (Table 1) and let $\mathrm{(}{\mu}^{\mathrm{\prime}}\mathrm{,}a\mathrm{,}b\mathrm{)}$ be the minimal doublemarked module in $\mathrm{M}$ onto which $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ is mapped by ${\mathrm{CG}}^{\mu \mathrm{\setminus}\mathrm{\{}a\mathrm{,}b\mathrm{\}}}$.
Reactions $a$ and $b$ are serial in $({\mu}^{\prime},a,b)$ if and only if they are serial in $(\mu ,a,b)$.
If reactions $a$ and $b$ are parallel in $({\mu}^{\prime},a,b)$, then they are also parallel in $(\mu ,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 $(\mu ,a,b)$, they may not be parallel in $({\mu}^{\prime},a,b)$. Figure 9 shows a counterexample illustrating this.
Corollary 4
Request a detailed protocolIf reactions $a$ and $b$ are strictly serial in $(\mu ,a,b)$, they are also strictly serial in $({\mu}^{\prime},a,b)$.
If reactions $a$ and $b$ are strictly parallel in $(\mu ,a,b)$, they are also strictly parallel in $({\mu}^{\prime},a,b)$.
If reactions $a$ and $b$ are serialparallel in $(\mu ,a,b)$, they are either strictly serial or serialparallel in $({\mu}^{\prime},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}_{\mu}$ is the same function that maps the rate constants $u$ and $v$ of these reactions in the minimal module ${\mu}^{\prime}$ onto the effective rate constant ${y}_{{\mu}^{\prime}}$. It then immediately follows from Equation 9 that the epistasis coefficient $\epsilon {y}_{\mu}$ between mutations affecting reactions $a$ and $b$ in the original module μ is the same as the epistasis coefficient $\epsilon {y}_{{\mu}^{\prime}}$ in the minimal module ${\mu}^{\prime}$. Now, if the reactions $a$ and $b$ are strictly parallel in $(\mu ,a,b)$, then, according to Corollary 4, these reactions are also strictly parallel in $({\mu}^{\prime},a,b)$. Hence, to demonstrate that $\epsilon {y}_{\mu}\le 0$ for any such doublemarked module $(\mu ,a,b)$, it is sufficient to show that $\epsilon {y}_{{\mu}^{\prime}}\le 0$ for all doublemarked modules $({\mu}^{\prime},a,b)$ that are minimal in $\mathcal{M}$ 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 doublemarked 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 nonpositive, irrespectively of the rate constants of any reactions, and similarly for the strictly serial reactions.
Generating topologies
Request a detailed protocolSince the number of topologically distinct minimal doublemarked 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 $f}_{2$), the idea is to identify the most connected minimal topologies (and the corresponding largest parametric families of $f}_{2$) 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 ${\mathcal{M}}^{\mathrm{b},\mathrm{io},\mathrm{IO}}$, ${\mathcal{M}}^{\mathrm{b},\mathrm{i},\mathrm{\varnothing}}$ and ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$ must be strictly parallel because the fully connected minimal topologies are strictly parallel (Figure 8). Similarly, all minimal modules in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{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 serialparallel. If the two reactions are serialparallel, 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 ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{\varnothing}}$, ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{I}}$, ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{\varnothing}}$, ${\mathcal{M}}^{\mathrm{i},\mathrm{i},\mathrm{I}}$, ${\mathcal{M}}^{\mathrm{i},\mathrm{i},\mathrm{\varnothing}}$ 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 serialparallel 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 $f}_{2$ which I then examine for the value of $\epsilon {y}_{\mu}$.
Proposition 5
Request a detailed protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ and $\mathrm{(}{\mu}^{\mathrm{\prime}}\mathrm{,}a\mathrm{,}b\mathrm{)}$ be two minimal doublemarked modules from the same topological class whose sets of reactions are ${R}_{\mu}$ and ${R}_{{\mu}^{\mathrm{\prime}}}$, respectively, and ${R}_{{\mu}^{\mathrm{\prime}}}\mathrm{=}{R}_{\mu}\mathrm{\setminus}\mathrm{\{}c\mathrm{\}}$ where $c\mathrm{\in}{R}_{\mu}\mathrm{\setminus}\mathrm{\{}a\mathrm{,}b\mathrm{\}}$.
If reactions $a$ and $b$ are strictly parallel in $(\mu ,a,b)$, they are also strictly parallel in $({\mu}^{\prime},a,b)$.
If reactions $a$ and $b$ are strictly serial in $(\mu ,a,b)$, they are also strictly serial in $({\mu}^{\prime},a,b)$.
Proof
Request a detailed protocolDenote the I/O metabolites in both modules μ and ${\mu}^{\prime}$ by 1 and 2. Since ${\mu}^{\prime}$ and μ are topologically identical except for ${\mu}^{\prime}$ lacking one reaction $c$, it must be true that ${\mathcal{P}}_{12}^{{\mu}^{\prime}}({a}_{1},{a}_{2},\mathrm{\dots},{\overline{b}}_{1},{\overline{b}}_{2},\mathrm{\dots})\subseteq {\mathcal{P}}_{12}^{\mu}({a}_{1},{a}_{2},\mathrm{\dots},{\overline{b}}_{1},{\overline{b}}_{2},\mathrm{\dots})$ for any reactions ${a}_{1},{a}_{2},\mathrm{\dots}$, ${b}_{1},{b}_{2},\mathrm{\dots}$ from ${R}_{{\mu}^{\prime}}$. In other words, there could only be fewer paths connecting the I/O metabolites within module ${\mu}^{\prime}$ 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 ${\mathcal{M}}^{\mathrm{b},\mathrm{io},\mathrm{IO}}$, ${\mathcal{M}}^{\mathrm{b},\mathrm{i},\mathrm{\varnothing}}$ and ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$, ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$) or if adding any reaction to it would make the marked reactions serialparallel.
Denote the sets of all doublemarked topologies minimal in class $\mathcal{M}$ where the marked reactions are strictly serial, strictly parallel and serialparallel by by ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{ser}}$, ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$ and ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$, respectively.
Definition 2
Request a detailed protocolTopology $\mathrm{(}R\mathrm{,}a\mathrm{,}b\mathrm{)}$ minimal in $\mathrm{M}$ is called a strictly serial generating topology in $\mathrm{M}$ if it is strictly serial (i.e. $\mathrm{(}R\mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{R}}_{\mathrm{M}}^{\mathrm{ser}}$) and either it is fully connected (i.e. $R\mathrm{=}{R}_{\mathrm{M}}$) or $\mathrm{(}R\mathrm{\cup}\mathrm{\{}c\mathrm{\}}\mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{R}}_{\mathrm{M}}^{\mathrm{sp}}$ for any reaction $c\mathrm{\in}{R}_{\mathrm{M}}\mathrm{\setminus}R$.
Definition 3
Request a detailed protocolTopology $\mathrm{(}R\mathrm{,}a\mathrm{,}b\mathrm{)}$ minimal in $\mathrm{M}$ is called a strictly parallel generating topology in $\mathrm{M}$ if it is strictly parallel (i.e. $\mathrm{(}R\mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{R}}_{\mathrm{M}}^{\mathrm{par}}$) and either it is fully connected (i.e. $R\mathrm{=}{R}_{\mathrm{M}}$) or $\mathrm{(}R\mathrm{\cup}\mathrm{\{}c\mathrm{\}}\mathrm{,}a\mathrm{,}b\mathrm{)}\mathrm{\in}{\mathrm{R}}_{\mathrm{M}}^{\mathrm{sp}}$ for any reaction $c\mathrm{\in}{R}_{\mathrm{M}}\mathrm{\setminus}R$.
Clearly, a topological class $\mathcal{M}$ 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 $\mathcal{M}$ by ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$ and I denote the set of all strictly parallel generating topologies for class $\mathcal{M}$ by ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{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 protocolIf $\mathrm{(}R\mathrm{,}a\mathrm{,}b\mathrm{)}$ is a strictly parallel topology minimal in the topological class $\mathrm{M}$, then there exists a strictly parallel generating topology $\mathrm{(}{R}_{g}\mathrm{,}a\mathrm{,}b\mathrm{)}$ in $\mathrm{M}$, such that $R\mathrm{\subseteq}{R}_{g}$. If $\mathrm{(}R\mathrm{,}a\mathrm{,}b\mathrm{)}$ is a strictly serial topology minimal in the topological class $\mathrm{M}$, then there exists a strictly serial generating topology $\mathrm{(}{R}_{g}\mathrm{,}a\mathrm{,}b\mathrm{)}$ in $\mathrm{M}$, such that $R\mathrm{\subseteq}{R}_{g}$.
Proof
Request a detailed protocolSuppose that $\mathcal{M}$ is one of the topological classes ${\mathcal{M}}^{\mathrm{b},\mathrm{io},\mathrm{IO}}$, ${\mathcal{M}}^{\mathrm{b},\mathrm{i},\mathrm{\varnothing}}$, or ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$. Since the fully connected minimal topology $({R}_{\mathcal{M}},a,b)$ is strictly parallel, it is a generating topology in $\mathcal{M}$. Then, according to Proposition 5, any topology $(R,a,b)$ minimal in $\mathcal{M}$ is strictly parallel, and Proposition 6 holds. By the same logic, Proposition 6 holds for the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$.
Suppose that $\mathcal{M}$ is one of the remaining topological classes ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{\varnothing}}$, ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{I}}$, ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{\varnothing}}$, ${\mathcal{M}}^{\mathrm{i},\mathrm{i},\mathrm{I}}$ or ${\mathcal{M}}^{\mathrm{i},\mathrm{i},\mathrm{\varnothing}}$. Then the fully connected minimal topology $({R}_{\mathcal{M}},a,b)$ is serialparallel. Suppose that $(R,a,b)$ is strictly parallel. Then $R$ must be a strict subset of ${R}_{\mathcal{M}}$, so that the set $C={R}_{\mathcal{M}}\setminus R$ is not empty. Then, either $(R,a,b)\in {\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$ or $(R,a,b)\notin {\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$. If $(R,a,b)\in {\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$, the Proposition 6 holds. Suppose that $(R,a,b)\notin {\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$. This implies that there exists a reaction ${c}_{1}\in C$, such that ${R}_{1}=R\cup \{{c}_{1}\}$ and $({R}_{1},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$ ($({R}_{1},a,b)$ cannot be in ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{ser}}$ due to Proposition 5). There are three possibilities.
${R}_{1}={R}_{\mathcal{M}}$.
${R}_{1}\subset {R}_{\mathcal{M}}$ and $({R}_{1},a,b)\in {\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$.
${R}_{1}\subset {R}_{\mathcal{M}}$ and $({R}_{1},a,b)\notin {\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$.
Option (a) is not possible since $({R}_{1},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$ while $({R}_{\mathcal{M}},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$. Option (b) implies that the Proposition 6 holds. Option (c) implies that there exists a reaction ${c}_{2}\in C\setminus \{{c}_{1}\}$, such that ${R}_{2}={R}_{1}\cup \{{c}_{2}\}$ and $({R}_{2},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$, and we have the same three possibilities for $R}_{2$ 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 $\mathcal{M}$ is straightforward because all minimal topologies within $\mathcal{M}$ can be produced by removing reactions from the unique fully connected topology minimal in $\mathcal{M}$ 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 fourletter 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, $\mathrm{i}\mathrm{o},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{1}$; $\mathrm{i}\mathrm{o},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{2}$, etc.
Topological relationship between reactions and epistasis
Request a detailed protocolEach strictly serial and strictly parallel generating topology in a given class $\mathcal{M}$ (listed in Table 2 and Table 3) is produced by removing reactions from the fully connected topology minimal in $\mathcal{M}$ (shown in Figure 8). This implies that the parametric family of function $f}_{2$ that corresponds to any generating topology is obtained from Equation 34 through Equation 35 by setting some parameters ${w}_{ij}$ to zero. In other words, these parametric families are known. Next, I prove Proposition 7 that shows that $\epsilon {y}_{\mu}\le 0$ for every member of every parameteric family of $f}_{2$ 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 $f}_{2$ that corresponds to any strictly parallel minimal topology belongs to the parametric family that corresponds to some strictly parallel generating topology; and any function $f}_{2$ 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 $\epsilon {y}_{\mu}\le 0$ for any minimal strictly parallel topology and that $\epsilon {y}_{\mu}\ge 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 $f}_{2$ that corresponds to some minimal strictly parallel module, and similarly for strictly serial modules.
Proposition 7
Request a detailed protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ be a minimal doublemarked module in the topological class $\mathrm{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 ${\delta}^{A}\mathit{}u$, and mutation B perturbs only reaction $b$ by ${\delta}^{B}\mathit{}v$, such that $\mathrm{\left}{\delta}^{A}\mathit{}u\mathrm{\right}\mathrm{\ll}\mathrm{1}$, $\mathrm{\left}{\delta}^{B}\mathit{}v\mathrm{\right}\mathrm{\ll}\mathrm{1}$. If reactions $a$ and $b$ are strictly parallel in $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$, then $\epsilon \mathit{}y\mathrm{\le}\mathrm{0}$.
Proof
Request a detailed protocolAccording to Proposition 6, the topology of module $(\mu ,a,b)$ can be produced by removing reactions from some strictly parallel generating topology $({R}_{g},a,b)$. Therefore, the function $f}_{2$ that maps $u$ and $v$ in this module onto its effective rate constant $y$ belongs to the parametric family that corresponds to $({R}_{g},a,b)$. According to Equation 55,
where
According to Theorem 1, $0\le {C}_{u}\le 1$ and $0\le {C}_{v}\le 1$. And since all $y>0$, $u>0$, $v>0$, to prove Proposition 7, it is sufficient to show that $\frac{{\partial}^{2}{f}_{2}}{\partial u\partial v}\le 0$ for any member of any parametric family that corresponds to generating topologies listed in Table 2.
Generating topologies $\mathrm{b}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{IO}\mathrm{,}\mathrm{F}$ and $\mathrm{b}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}\mathrm{F}$ (Figure 8). According to Equation 34 and Equation 35, $y={f}_{2}(u,v)=u+\varphi (v)$, which implies that $\epsilon y=0$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{IO}\mathrm{,}\mathrm{F}$ (Figure 8). According to Equation 37,
where $D=Euv+Fu+Gv+B$, $A={w}_{42}/{K}_{13}+{w}_{32}/{K}_{14}$, $B={w}_{32}{w}_{42}+{w}_{32}{w}_{43}+{w}_{34}{w}_{42}$, $E=1/\left({K}_{13}{K}_{14}\right)$, $F=\left({w}_{42}+{w}_{43}\right)/{K}_{13}$, $G=\left({w}_{32}+{w}_{34}\right)/{K}_{14}$. Therefore,
Generating topology $\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{\varnothing}\mathrm{,}\mathrm{P}$ (Figure 10). According to Equation 38, $y={f}_{2}(u,v)={w}_{12}+\varphi (u)+\psi (v)$, which implies $\epsilon y=0$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}\mathrm{P}$ (Figure 11). Notice that metabolite 4 together with reactions $1\leftrightarrow 4$, $a$ and $b$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 3 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{b},\mathrm{io},\mathrm{IO}}$. Denote the effective reaction rate of module ${\mu}^{\prime}$ by ${y}^{\prime}$. Therefore, $\epsilon {y}^{\prime}=0$, as shown above. Since module ${\mu}^{\prime}$ is contained in μ, by Theorem 1, $\epsilon y\le 0$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{1}}$ (Figure 12). According to Property 1 of the CGP (Box 1), module μ can be coarsegrained by first eliminating metabolite 3. In the resulting module ${\mu}^{\prime}$, mutation A perturbs only the rate constant ${u}^{\prime}$ of the effective reaction ${a}^{\prime}\equiv 1\leftrightarrow 2$ (by Properties 2 and 4 of the CGP). Then, according to Equation 50 and Theorem 1, $\left{\delta}^{A}{u}^{\prime}\right\ll 1$. The doublemarked module $({\mu}^{\prime},{a}^{\prime},b)$ is minimal in the topological class ${\mathcal{M}}^{\mathrm{b},\mathrm{i},\mathrm{\varnothing}}$ which implies that $\epsilon y=0$, as shown above.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{2}}$ (Figure 12). Module μ can be coarsegrained by first eliminating metabolite 5, which will result in a doublemarked module $({\mu}^{\prime},a,{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{1}$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{3}}$ (Figure 12). Notice that metabolites 4 and 5 together with reactions $a$, $b$, $1\leftrightarrow 4$, $1\leftrightarrow 5$, $3\leftrightarrow 4$ and $3\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 3 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{b},\mathrm{i},\mathrm{\varnothing}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},\mathrm{P}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}{\mathrm{P}}_{\mathrm{1}}$ (Figure 13). Notice that metabolites 4 and 5 together with reactions $a$, $b$, $1\leftrightarrow 3$, $1\leftrightarrow 4$, $1\leftrightarrow 5$ and $4\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 3 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},\mathrm{P}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}{\mathrm{P}}_{\mathrm{2}}$ (Figure 13). Notice that metabolite 5 together with reactions $a$, $b$, and $4\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 3 and 4 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{b},\mathrm{io},\mathrm{IO}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},\mathrm{P}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{1}}$ (Figure 14). According to Equation 35, $y={f}_{2}(u,v)={x}_{12}+\varphi (u)+\psi (v)$, which implies $\epsilon y=0$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{2}}$ (Figure 14). Notice that metabolites 5 and 6 together with reactions $a$, $b$, $3\leftrightarrow 5$, $3\leftrightarrow 6$, $4\leftrightarrow 5$ and $5\leftrightarrow 6$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 3 and 4 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{b},\mathrm{i},\mathrm{\varnothing}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},\mathrm{P}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{3}}$ (Figure 14). Module μ can be coarsegrained by first eliminating metabolites 4 and 6, which will result in a doublemarked module $({\mu}^{\prime},{a}^{\prime},{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{IO}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{1}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{4}}$ (Figure 14). Module μ can be coarsegrained by first eliminating metabolite 4, which will result in a doublemarked module $({\mu}^{\prime},{a}^{\prime},b)$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{\varnothing}}$ with a strictly parallel generating topology $\mathrm{io},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{3}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{1}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{5}}$ (Figure 14). Module μ can be coarsegrained by first eliminating metabolite 6, which will result in a doublemarked module $({\mu}^{\prime},a,{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{i},\mathrm{i},\mathrm{I}}$ with a strictly parallel generating topology $\mathrm{i},\mathrm{i},\mathrm{I},{\mathrm{P}}_{1}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{1}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{6}}$ (Figure 14). Notice that metabolites 4 and 6 together with reactions $a$, $b$, $3\leftrightarrow 5$, $3\leftrightarrow 6$, $4\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 3 and 5 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{\varnothing}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},\mathrm{P}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{P}}_{\mathrm{7}}$ (Figure 14). Using Equation 42, I show in Appendix 4 that
where ${A}_{v}$, ${B}_{v}$, ${E}_{u}$, ${F}_{u}$ are all nonnegative constants and $\beta \le 0$.
Proposition 8
Request a detailed protocolLet $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$ be a minimal doublemarked module in the topological class $\mathrm{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 ${\delta}^{A}\mathit{}u$, and mutation B perturbs only reaction $b$ by ${\delta}^{B}\mathit{}v$, such that $\mathrm{\left}{\delta}^{A}\mathit{}u\mathrm{\right}\mathrm{\ll}\mathrm{1}$, $\mathrm{\left}{\delta}^{B}\mathit{}v\mathrm{\right}\mathrm{\ll}\mathrm{1}$. If reactions $a$ and $b$ are strictly serial in $\mathrm{(}\mu \mathrm{,}a\mathrm{,}b\mathrm{)}$, then $\epsilon \mathit{}y\mathrm{\ge}\mathrm{1}$.
Proof
Request a detailed protocolThe logic of the proof is the same as for Proposition 7, that is, it is sufficient to show that $\epsilon y\ge 1$ for any doublemarked module $(\mu ,a,b)$ with any strictly serial generating topology listed in Table 3.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{I}\mathrm{,}\mathrm{F}$ (Figure 8). According to Equation 36,
where $D=u/{K}_{13}+v$. Therefore,
Substituting these expressions into Equation 68, I obtain
because $y\ge uv/D$ according to Equation 73.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{io}\mathrm{,}\mathrm{\varnothing}\mathrm{,}\mathrm{S}$ (Figure 10). According to Property 1 of the CGP (Box 1), module μ can be coarsegrained by first eliminating metabolite 3. In the resulting module ${\mu}^{\prime}$, mutation A perturbs only the rate constant ${u}^{\prime}$ of the effective reaction ${a}^{\prime}\equiv 1\leftrightarrow 4$ (by Properties 2 and 4 of the CGP). Then, according to Equation 50 and Theorem 1, $\left{\delta}^{A}{u}^{\prime}\right\ll 1$. The doublemarked module $({\mu}^{\prime},{a}^{\prime},b)$ is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$ which implies that $\epsilon y\ge 1$, as shown above.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}{\mathrm{S}}_{\mathrm{1}}$ (Figure 11). Notice that metabolite 3 together with reactions $a$, $b$, and $1\leftrightarrow 4$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 4 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$. Therefore, if the effective reaction rate of module ${\mu}^{\prime}$ is ${y}^{\prime}$, $\epsilon {y}^{\prime}\ge 1$, as shown above. According to Equation 50, Equation 51 and Theorem 1, $\left{\delta}^{A}{y}^{\prime}\right\ll 1$, $\left{\delta}^{B}{y}^{\prime}\right\ll 1$. Since module ${\mu}^{\prime}$ is contained in μ, by Theorem 1, $\epsilon y\ge 1$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}{\mathrm{S}}_{\mathrm{2}}$ (Figure 11). Module μ can be coarsegrained by first eliminating metabolite 4, which results in a doublemarked module $({\mu}^{\prime},a,{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{io},\mathrm{\varnothing},\mathrm{S}$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{S}}_{\mathrm{1}}$ (Figure 12). Module μ can be coarsegrained by first eliminating metabolites 4 and 5, which results in a doublemarked module $({\mu}^{\prime},a,{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{io},\mathrm{\varnothing},\mathrm{S}$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{S}}_{\mathrm{2}}$ (Figure 12). Notice that metabolites 3 and 4 together with reactions $a$, $b$, $1\leftrightarrow 4$, $1\leftrightarrow 5$ and $3\leftrightarrow 4$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 5 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{\varnothing}}$ with the strictly serial generating topology $\mathrm{io},\mathrm{io},\mathrm{\varnothing},\mathrm{S}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},{\mathrm{S}}_{1}$.
Generating topology $\mathrm{io}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{S}}_{\mathrm{3}}$ (Figure 12). Notice that metabolites 3 and 5 together with reactions $a$, $b$, $1\leftrightarrow 4$, $3\leftrightarrow 4$, and $3\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 4 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{\varnothing}}$ with the strictly serial generating topology $\mathrm{io},\mathrm{io},\mathrm{\varnothing},\mathrm{S}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},{\mathrm{S}}_{1}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}{\mathrm{S}}_{\mathrm{1}}$ (Figure 13). Notice that metabolite 3 together with reactions $a$, $b$, and $4\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 4 and 5 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{io},\mathrm{I}}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},{\mathrm{S}}_{1}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{I}\mathrm{,}{\mathrm{S}}_{\mathrm{2}}$ (Figure 13). Notice that metabolites 3 and 5 together with reactions $a$, $b$, $1\leftrightarrow 3$, $1\leftrightarrow 4$, and $1\leftrightarrow 5$ form a doublemarked module $({\mu}^{\prime},a,b)$ whose I/O metabolites are 1 and 4 and which is minimal in the topological calss ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{I}}$ with the strictly serial generating topology $\mathrm{io},\mathrm{i},\mathrm{I},{\mathrm{S}}_{2}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{i},\mathrm{I},{\mathrm{S}}_{1}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{S}}_{\mathrm{1}}$ (Figure 14). Module μ can be coarsegrained by first eliminating metabolites 5 and 6, which results in a doublemarked module $({\mu}^{\prime},a,{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{io},\mathrm{i},\mathrm{I}}$ with the strictly serial generating topology $\mathrm{io},\mathrm{i},\mathrm{I},{\mathrm{S}}_{1}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{io},\mathrm{\varnothing},\mathrm{S}$.
Generating topology $\mathrm{i}\mathrm{,}\mathrm{i}\mathrm{,}\mathrm{\varnothing}\mathrm{,}{\mathrm{S}}_{\mathrm{2}}$ (Figure 14). Module μ can be coarsegrained by first eliminating metabolite 6, which results in a doublemarked module $({\mu}^{\prime},a,{b}^{\prime})$ that is minimal in the topological class ${\mathcal{M}}^{\mathrm{i},\mathrm{i},\mathrm{I}}$ with the strictly serial generating topology $\mathrm{i},\mathrm{i},\mathrm{I},{\mathrm{S}}_{1}$. The rest of the proof for this topology is analogous to that for the topology $\mathrm{io},\mathrm{io},\mathrm{\varnothing},\mathrm{S}$.
Proof of Theorem 2
Request a detailed protocolAccording to Proposition 3, the coarsegraining procedure ${\mathrm{CG}}^{\mu \setminus \{a,b\}}$ maps the doublemarked module $(\mu ,a,b)$ onto a doublemarked module $({\mu}^{\prime},a,b)$ that is minimal in the same topological class as $(\mu ,a,b)$, and the rates $u$, $v$ of reactions $a$, $b$ in ${\mu}^{\prime}$ are given by linear relations in Equation 32 and Equation 33. Clearly, $\left{\delta}^{A}u\right\ll 1$ and $\left{\delta}^{B}v\right\ll 1$. Furthermore, none of the other reaction rates ${w}_{ij}$ in ${\mu}^{\prime}$ depend on ξ or η, so that ${\delta}^{A}{w}_{ij}=0$ and ${\delta}^{B}{w}_{ij}=0$ for all ${w}_{ij}$ other than $u$ and $v$, and $\epsilon {w}_{ij}=0$ for all ${w}_{ij}$ including $u$ and $v$. It then follows from Proposition 3 that $\epsilon {y}_{\mu}=\epsilon {y}_{{\mu}^{\prime}}$.
Now, according to Corollary 4, if reactions $a$ and $b$ are strictly parallel in $(\mu ,a,b)$, they are also strictly parallel in $({\mu}^{\prime},a,b)$. Therefore, by Proposition 7, $\epsilon {y}_{{\mu}^{\prime}}\le 0$. Analogously, if reactions $a$ and $b$ are strictly serial in $(\mu ,a,b)$, they are also strictly serial in $({\mu}^{\prime},a,b)$. Therefore, by Proposition 8, $\epsilon {y}_{{\mu}^{\prime}}\ge 1$.
Sensitivity of Theorem 1 and Theorem 2 with respect to the magnitude of mutational effects
Request a detailed protocolAccording to Proposition 2, function $f}_{1$ 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 ${\mathcal{M}}^{\mathrm{b}}$, function f_{1} 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 ${\mathcal{M}}^{\mathrm{io}}$ and ${\mathcal{M}}^{\mathrm{i}}$, I generated 1000 minimal singlemarked 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}_{\mu}$ onto the higherlevel phenotype ${y}_{\nu}$ via the same function f_{1} (see Proposition 2).
To this end, I drew each ${x}_{ij}$ ($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 ${K}_{ij}$ ($i<j$) as a ratio of two random numbers from an exponential distribution with mean 1. As a result, the distribution of nonzero ${x}_{ij}$ values had the interdecile range of $(5.7\times {10}^{2},3.91)$ with median 0.65.
I denote the effective rate constant of the reaction that represents the lowerlevel module μ by $\xi \equiv {y}_{\mu}$. In modules from the topological class ${\mathcal{M}}^{\mathrm{io}}$, it is reaction $1\leftrightarrow 3$ and in modules from the topological class ${\mathcal{M}}^{\mathrm{i}}$, it is the reaction $3\leftrightarrow 4$. I perturbed ξ by two mutations A and B with relative effects ${\delta}^{A}\xi $ and ${\delta}^{B}\xi $ and epistasis $\epsilon \xi $. I chose nine different pairs of mutational effects $({\delta}^{A}\xi ,{\delta}^{B}\xi ):(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 $\epsilon \xi $ ranging from −1 to 2 with an increment of 0.2. Since the rate constant ${\xi}^{AB}$ of the double mutant cannot be negative, I skipped those combinations of perturbations and epistasis values for which ${\delta}^{A}\xi +{\delta}^{B}\xi +2\left(\epsilon \xi \right)\left({\delta}^{A}\xi \right)\left({\delta}^{B}\xi \right)<1$. I then computed the resulting values ${\delta}^{A}{y}_{\nu}$, ${\delta}^{B}{y}_{\nu}$ and $\epsilon {y}_{\nu}$ at the level of the effective rate constant ${y}_{\nu}$ of the higherlevel module ν.
Using these data, I inferred the function $\varphi $ that maps lowerlevel epistasis $\epsilon \xi $ onto higherlevel epistasis $\epsilon {y}_{\nu}$, as follows. For any minimal singlemarked module from the topological classes ${\mathcal{M}}^{\mathrm{io}}$ or ${\mathcal{M}}^{\mathrm{i}}$, the effective rate constant ${y}_{\nu}$ can be written as
where $D=\stackrel{~}{C}\xi +\stackrel{~}{D}$ and $\stackrel{~}{A}={x}_{32}$, $\stackrel{~}{B}=0$, $\stackrel{~}{C}=1/{K}_{13}$ , $\stackrel{~}{D}={x}_{32}$ for modules from the topological class ${\mathcal{M}}^{\mathrm{io}}$ (see Equation 30), and $\stackrel{~}{A}=\left({x}_{13}+{x}_{14}\right)\left({x}_{42}+{x}_{32}/{K}_{34}\right)$, $\stackrel{~}{B}=({x}_{31}+{x}_{32}){x}_{14}{x}_{42}+({x}_{41}+{x}_{42}){x}_{13}{w}_{32}$, $\stackrel{~}{C}=({x}_{31}+{x}_{32})/{K}_{34}+({x}_{41}+{x}_{42})$ , $\stackrel{~}{D}=({x}_{31}+{x}_{32})({x}_{41}+{x}_{42})$ for modules from the topological class ${\mathcal{M}}^{\mathrm{i}}$ (see Equation 31). Therefore, for any perturbation $\delta \xi $, we have
Since ${\delta}^{AB}\xi $ is a linear function of $\epsilon \xi $, ${\delta}^{AB}{y}_{\nu}$ is a hyperbolic function of $\epsilon \xi $. Therefore, $\epsilon {y}_{\nu}$ is also a hyperbolic function of $\epsilon \xi $,
where constants $a$, $b$ and $c$ depend on the parameters of module ν and on the mutational effect sizes ${\delta}^{A}\xi $ and ${\delta}^{B}\xi $. 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 $\varphi $ has a fixed point $\overline{\epsilon}$, 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 $\varphi $ has a fixed point $\overline{\epsilon}$ if the discriminant $d={(ac)}^{2}4(bac)$ is positive. In this case, I designate $\overline{\epsilon}$ as the one of two roots $1/2\left(ac\pm \sqrt{d}\right)$ 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 $\varphi $ at $\overline{\epsilon}$ with 1.
According to Proposition 6, function $f}_{2$ 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 $f}_{2$ 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 10^{4} doublemarked modules $(\mu ,a,b)$ with each of the strictly serial and strictly parallel topologies with random parameters. I drew ${x}_{ij}$ and ${K}_{ij}$ as described above. I chose the same nine pairs of mutational effects $({\delta}^{A}\xi ,{\delta}^{B}\eta )$ 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 ${\delta}^{A}{y}_{\mu}$ and/or ${\delta}^{B}{y}_{\mu}$ at the level of the whole module were too small, which resulted in numerical instabilities. To avoid them, I calculated epistasis $\epsilon {y}_{\mu}$ only for cases where the effects of both mutations ${\delta}^{A}{y}_{\mu}$ and ${\delta}^{B}{y}_{\mu}$ exceeded the precision threshold of ${10}^{5}$. As a result, I evaluated epistasis in less than 10^{4} 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 $\epsilon {y}_{\mu}$ 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 (MalikSheriff 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 protocolI simplified and modified the model by (a) fixing the concentrations of ATP, ADP, AMP, NADPH, NADP, NADH, NAD at their steadystate values given in Table V of Chassagnole et al., 2002 and (b) removing dilution by growth. I then created four models of submodules 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 submodule 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 steadystate values given in Table 4. I defined the flux through the submodule as the flux toward the output metabolite contributed by the submodule (Table 4). This flux is the equivalent of the quantitative phenotype ${y}_{\mu}$ of a module in the analytical model. In addition, I made the following modifications specific to individual submodules.
In the FULL model, the stoichiometry of the PTS reaction was changed to
$[\text{Ext glu}]+[\text{pep}]\leftrightarrow [\text{g6p}]+[\text{pyr}]$and the value of the constant ${K}_{\mathrm{PTS},\mathrm{a1}}$ was set to 0.02 mM, based on the values found in the literature (Stock et al., 1982; Natarajan and Srienc, 1999).
In all models other than FULL, the extracellular compartment was deleted.
In all models, the concentrations of the I/O metabolites were set to values shown in Table 4, which are the steadystate 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 protocolI calculate the first and secondorder flux control coefficients (FCC) ${C}_{i}$ and ${H}_{ij}$ for flux $J$ with respect to reactions $i$ and $j$ as follows (see Equation 45 and Equation 46). I perturb the ${r}_{\mathrm{max},i}$ of reaction $i$ by factor between 0.75 and 1.25 (10 values in a uniformlyspaced grid), such that $\delta {r}_{\mathrm{max},i}\in [0.25,0.25]$. Then, I obtain the steadystate flux ${J}^{\prime}$ in each perturbed model and calculate the flux perturbations $\delta J={J}^{\prime}/{J}^{0}1$, where ${J}^{0}$ is the corresponding flux in the unperturbed model. Then, to obtain ${C}_{i}$ and ${H}_{ii}$, I fit the linear model
by least squares. If the estimated value of ${C}_{i}$ was below ${10}^{4}$ for a given submodule, I set ${C}_{i}$ to zero and exclude this reaction from further consideration in that submodule because it does not affect flux to the degree that is accurately measurable. If the estimated value of ${H}_{ii}$ is below ${10}^{4}$, I set ${H}_{ii}$ to zero.
To calculate the nondiagonal secondorder control coefficients ${H}_{ij}$, I create a $4\times 4$ grid of perturbations of $\delta {r}_{\mathrm{max},i}$ and $\delta {r}_{\mathrm{max},j}$ and calculate the resulting flux perturbations $\delta J$ (16 perturbations total). Since ${C}_{i}$, ${C}_{j}$, ${H}_{ii}$ and ${H}_{jj}$ are known, I obtain ${H}_{ij}$, by regressing
against
If the estimated value of ${H}_{ij}$ is below ${10}^{4}$, I set ${H}_{ij}$ to zero. I estimate the epistasis coefficient $\epsilon J$ between mutations affecting reactions $i$ and $j$ as
Establishing the topological relationships between pairs of reactions
Request a detailed protocolTo establish the topological relationship (strictly serial, strictly parallel, or serialparallel) 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 serialparallel.
Appendix 1
Proof of Equation 24
First, the terms in Equation 24 can be rearranged as follows
Next, I will demonstrate the validity of Equation 75 by induction. It is clear that for $E=\{k\}\subset {A}_{\mu}$, Equation 24 reduces to Equation 12. Now suppose that ${n}_{E}>1$ and that there exists a metabolite $k\in E$, such that Equation 75 holds for the subset ${E}^{\prime}=E\setminus \{k\}$, that is,
I will now show that then Equation 75 also holds for $E$. To do so, I use the definition of ${x}_{ij}^{E}$ (Equation 19) and apply Equation 76 to expand terms ${x}_{ij}^{{E}^{\prime}}$ and ${x}_{ik}^{{E}^{\prime}}$ as follows.
Recalling that $E={E}^{\prime}\setminus \{k\}$, I rewrite Equation 75 as
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 $\mathrm{\ell}\in {E}^{\prime}$,
To show that Equation 79 holds, I first use Equation 19 and Equation 21 to express ${x}_{kj}^{E\setminus \{k\}}$ and ${D}_{k}^{E\setminus \{k\}}$ in terms of effective reaction rates after the elimination of the metabolite set $E\setminus \{k,\mathrm{\ell}\}$,
which imply that
Finally, it follows from Equation 80 that ${D}_{\mathrm{\ell}}^{E\setminus \{k,\mathrm{\ell}\}}{D}_{k}^{E\setminus \{k\}}={D}_{k}^{E\setminus \{k,\mathrm{\ell}\}}{D}_{\mathrm{\ell}}^{E\setminus \{\mathrm{\ell}\}},$ which completes the proof of Equation 79. Thus, Equation 75 and equivalently Equation 24 hold for any metabolite subset $E\subseteq {A}_{\mu}$.
Appendix 2
Existence of a simple path that contains a given reaction
Lemma 1
Let μ be a module with the reaction set ${R}_{\mu}$. Then, for any reaction $a\mathrm{\in}{R}_{\mu}$, there exists a simple path ${p}_{\mathrm{12}}\mathit{}\mathrm{(}a\mathrm{)}$ within μ that connects the I/O metabolites and contains reaction $a$, that is, ${\mathrm{P}}_{\mathrm{12}}^{\mu}\mathit{}\mathrm{(}a\mathrm{)}\mathrm{\ne}\mathrm{\varnothing}$.
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=1\leftrightarrow j$. Since μ is a module, there exists a simple path $j\leftrightarrow {j}_{1}\leftrightarrow \mathrm{\cdots}\leftrightarrow 2$ that connects the internal metabolite j to the I/O metabolite 2. Therefore, the path $1\leftrightarrow j\leftrightarrow {j}_{1}\leftrightarrow \mathrm{\cdots}\leftrightarrow 2$ connects the I/O metabolites and contains reaction $a$.
Suppose that $a=i\leftrightarrow j$ is an internal reaction. To prove the statement, it is sufficient to show that there exists a pair of nonintersecting paths ${p}_{1i}^{\prime}$ and ${p}_{2i}^{\prime}$, such that one of them contains $a$ and the other does not. Since μ is a module, there exists a pair of nonintersecting paths ${p}_{1i}$ and ${p}_{2i}$ and a pair of nonintersecting paths ${p}_{1j}$ and ${p}_{2j}$ within module μ (I omitted superindex μ to simplify notations). There are two mutually exclusive possibilities. (i) Metabolite $j$ is contained in either of the paths ${p}_{1i}$ or ${p}_{2i}$ and/or metabolite $i$ is contained in either of the paths ${p}_{1j}$ or ${p}_{2j}$. (ii) Metabolite $j$ is not contained in either of the paths ${p}_{1i}$ or ${p}_{2i}$ and metabolite $i$ is not contained in either of the paths ${p}_{1j}$ or ${p}_{2j}$, that is, $j\notin {p}_{ui}$ and $i\notin {p}_{uj}$, $u=1,2$. It is trivial to construct the neccessary paths ${p}_{1i}^{\prime}$ and ${p}_{2i}^{\prime}$ in case (i).
Consider case (ii). If paths ${p}_{2j}$ and ${p}_{1i}$ do not intersect, then let ${p}_{1i}^{\prime}={p}_{1i}$ and
and the statement is true. Suppose paths ${p}_{2j}$ and ${p}_{1i}$ intersect. Then, among all metabolites that belong to both ${p}_{1i}$ and ${p}_{2j}$, let metabolite $k$ be the one closest to $j$ along the path ${p}_{2j}$ (see Figure). Then the segment ${p}_{kj}$ of path ${p}_{2j}$ and the path ${p}_{1i}$ do not intersect. Let
If ${p}_{1i}^{\prime \prime}$ and ${p}_{2i}$ do not intersect, then the lemma is true. If ${p}_{1i}^{\prime \prime}$ and ${p}_{2i}$ do intersect, this intersetion can only occur within the segment ${p}_{kj}$ of path ${p}_{1i}^{\prime \prime}$, excluding metabolites $k$ and $j$ (see Figure). This is because the remaining segment ${p}_{1k}$ of path ${p}_{1i}^{\prime \prime}$ is also a segment of ${p}_{1i}$, which, by assumption, does not intersect ${p}_{2i}$. Suppose that among all metabolites that belong to both the segment ${p}_{kj}$ of path ${p}_{1i}^{\prime \prime}$ and the path ${p}_{2i}$ metabolite $\mathrm{\ell}$ is the one closest to $j$ along the path ${p}_{1i}^{\prime \prime}$. Then let
The path ${p}_{2i}^{\prime}$ does not intersect the path ${p}_{1i}$ because its first segment ${p}_{2\mathrm{\ell}}$ belongs to path ${p}_{2i}$ and its second segment ${p}_{\mathrm{\ell}j}$ belongs to the segment ${p}_{kj}$ of path ${p}_{1i}^{\prime \prime}$ (and, as mentioned above, segment ${p}_{kj}$ does not intersect ${p}_{1i}$). 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)\in {\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$, that is, $(R,a,b)$ is a strictly serial generating topology in the topological class $\mathcal{M}$. Since $R\subseteq {R}_{\mathcal{M}}$ where ${R}_{\mathcal{M}}$ is the complete reaction set ${R}_{\mathcal{M}}$ for class $\mathcal{M}$, R can be discovered by sequentially removing reactions from ${R}_{\mathcal{M}}$. The same logic holds for strictly parallel generating topologies. The following algorithm implements this idea.
Define function generate_topology_list. This function takes a topology $(R,a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$ as input and returns a new list of topologies $L$ as output, which is produced as follows. Initialize $L=\mathrm{\varnothing}$. For every reaction ${c}_{i}\in R\setminus \{a,b\}$, construct the reaction subset ${R}_{i}=R\setminus \{{c}_{i}\}$ and use Definition 1 to test whether ${R}_{i}$ corresponds to a valid module. If ${R}_{i}$ corresponds to a module, add $({R}_{i},a,b)$ to list $L$; otherwiese, discard. It can be proven that, as long as $(R,a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$, there exists at least one ${c}_{i}\in R$, such that ${R}_{i}$ corresponds to a module, that is, $L\ne \mathrm{\varnothing}$. Return list $L$.
Initialization.
Pick a topological class $\mathcal{M}$.
Test whether $({R}_{\mathcal{M}},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{ser}}$. If so, ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}=\left\{({R}_{\mathcal{M}},a,b)\right\}$ and ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}=\mathrm{\varnothing}$. Return ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$, ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$.
Test whether $({R}_{\mathcal{M}},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$. If so, ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}=\left\{({R}_{\mathcal{M}},a,b)\right\}$ and ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}=\mathrm{\varnothing}$. Return ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$, ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$.
Set ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}=\mathrm{\varnothing}$, ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}=\mathrm{\varnothing}$. Use function generate_topology_list with $({R}_{\mathcal{M}},a,b)$ as input and obtain the list of reaction sets $L$. Proceed to Step 3 with list $L$.
Take list $L=(({R}_{1},a,b),({R}_{2},a,b),\mathrm{\dots},({R}_{k},a,b))$ as input. Set ${L}^{\prime}=\mathrm{\varnothing}$. Proceed to Step a with $i=1$.
Test whether $({R}_{i},a,b)$ belongs to ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$, ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{ser}}$ or ${\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$. Choose one of the alternatives b, c, or d.
$({R}_{i},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{par}}$. If $({R}_{i},a,b)$ is not in the set ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$ and if $({R}_{i}\cup \{c\},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$ for all $c\in {R}_{\mathcal{M}}\setminus {R}_{i}$, then add $({R}_{i},a,b)$ to ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$. Proceed to Step e.
$({R}_{i},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{ser}}$. If $({R}_{i},a,b)$ is not in the set ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$ and if $({R}_{i}\cup \{c\},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$ for all $c\in {R}_{\mathcal{M}}\setminus {R}_{i}$, then add $({R}_{i},a,b)$ to ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$. Proceed to Step e.
$({R}_{i},a,b)\in {\mathcal{R}}_{\mathcal{M}}^{\mathrm{sp}}$. Use function generate_topology_list with $({R}_{i},a,b)$ as input and obtain the list of reaction sets ${L}_{i}$. Replace ${L}^{\prime}$ with ${L}^{\prime}\cup {L}_{i}$. Proceed to Step e.
If $i=k$, proceed to Step f. Otherwise, proceed to Step a with $i+1$.
If ${L}^{\prime}\ne \mathrm{\varnothing}$, then proceed to Step 3 with ${L}^{\prime}$ as input. Otherwise, return ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{ser}}$, ${\mathcal{G}}_{\mathcal{M}}^{\mathrm{par}}$.
Appendix 4
Derivation of Equation 72
In this case, Equation (42) simplify to
where
Notice that the only effective rate constant that depends on $u$ is $W}_{34$, and $\frac{\partial {W}_{34}}{\partial u}=1$. Thus, it is easy to differentiate $y$, given by Equation 35, with respect to $u$, if we isolate the term $W}_{34$ in both the numerator and the denomintor,
where
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
where
The only effective activity that depends on $v$ is $V}_{56$, and $\frac{\mathrm{\partial}{V}_{56}}{\mathrm{\partial}v}=1$. Thus, isolating the term $V}_{56$ (and recalling that $V}_{65}={V}_{56}/{K}_{56$), we obtain the following expression for $y$, which is easy to differentiate with respect to $v$.
where
Using symbolic computation it is possible to show that (see Mathematica notebook [Supplementary file 1] and Mathematica notebook pdf [Supplementary file 2], p. 2)
Also notice that, any module with the reaction set $\mathrm{i},\mathrm{i},\mathrm{\varnothing},{\mathrm{P}}_{7}$ 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
Analogously, differentiating Equation 82 with respect to $v$, I obtain
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
where all coefficients
and
are independent of $u$ and $v$ and are nonnegative. Similarly (see Mathematica notebook (Supplementary file 1) and Mathematica notebook pdf (Supplementary file 2), p. 4),
where
We can now obtain the second derivative $\frac{{\partial}^{2}y}{\partial u\partial v}$, taking into account Equation 86. Alternatively, we can obtain $\frac{{\partial}^{2}y}{\partial u\partial v}$ by differentiating $\frac{\partial y}{\partial v}$ 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,
where β is independent of $u$ and $v$. Thus, according to Equation 84, Equation 86,
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,
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 wholecell 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 steadystate concentrations of internal metabolites (see e.g. Savageau, 1976). The steadystate 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 steadystate fluxes are subject only to the massbalance 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 tradesoff 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 firstorder kinetics described by my model are a special case. So, let us apply FBA to a metabolic module $\mu =({A}_{\mu},{\overrightarrow{x}}_{\mu})$ 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 ${A}_{\mathrm{in}}\subset {A}_{\mu}$ are adjacent to the external nutrient and metabolites in the set ${A}_{\mathrm{out}}\subset {A}_{\mu}$ are adjacent to the biomass. Then reactions $1\leftrightarrow i$ for $i\in {A}_{\mathrm{in}}$, that is, those that convert the nutrient into intermediate metabolites, are the ‘uptake reactions’. And the reactions $i\leftrightarrow 2$ for $i\in {A}_{\mathrm{out}}$, 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 ${v}_{ij}$ and let ${J}_{\mathrm{in}}={\sum}_{i\in {A}_{\mathrm{in}}}{v}_{1i}$ and ${J}_{\mathrm{out}}={\sum}_{i\in {A}_{\mathrm{out}}}{v}_{i2}$ be the input and output fluxes, respectively. To study epistasis within the FBA framework, we would need to obtain ${J}_{\mathrm{out}}$.
The assumption of my model that all reaction kinetics are firstorder translates into the FBA formalism as the fact that the elements of the stoichiometry matrix $\overrightarrow{S}$ take values −1, +1 or 0 and there are no capacity constraints on the rates of reactions. Then, the massbalance equation $\overrightarrow{S}\phantom{\rule{thinmathspace}{0ex}}\overrightarrow{v}=\overrightarrow{0}$ leads to the equality ${J}_{\mathrm{in}}={J}_{\mathrm{out}}$, 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 steadystate flux $J$ through the module depends in general on the internal structure of the module (through the effective rate constant ${y}_{\mu}$) and is given by $J={y}_{\mu}\left({S}_{1}{S}_{2}/{K}_{12}\right)$ for any concentrations $S}_{1$ and $S}_{2$ of the external nutrient and the biomass. However, in the degenerate case where all uptake reactions are irreversible (i.e., ${K}_{1i}=\mathrm{\infty}$ for all $i\in {A}_{\mathrm{in}}$), $J={J}_{\mathrm{in}}={\sum}_{i\in {A}_{\mathrm{in}}}{x}_{1i}{S}_{1}$, 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 ${S}_{i}$ 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 ${A}_{i}={S}_{i}vN$, $i=1,\mathrm{\dots},n1$.
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 $i\leftrightarrow j$, $i<j$ be ${\varphi}_{ij}$ (with ${\sum}_{i=1}^{n1}{\sum}_{j>i}{\varphi}_{ij}=1$), so that its concentration inside cells is ${\varphi}_{ij}p$. If the specific forward and reverse activities of this enzyme are ${a}_{ij}$ and ${a}_{ji}$ (so that ${a}_{ij}$ obey the Haldane equalities and ${a}_{ii}=0$), then the total forward and reverse activities are ${x}_{ij}={a}_{ij}{\varphi}_{ij}p$ and ${x}_{ji}={a}_{ji}{\varphi}_{ij}p$.
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 $S}_{1$ stays constant. In other words, I model early exponential growth. Then the dynamics of metabolite and protein abundances are governed by equations
where ${D}_{i}={\sum}_{i=1}^{n}{x}_{ij}$.
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
where $\lambda ={\sum}_{i=2}^{n1}{x}_{in}{S}_{i}$ is the steadystate 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 $\lambda {S}_{i}\approx 0$. Then Equation 90 defines the same steady state for module μ as described in the section Network coarsegraining. Therefore, module μ can be replaced with the effective activity ${x}_{\mu}$ between nutrients and biomass, yielding $\lambda ={x}_{\mu}$.
Data availability
All data generated or analyzed during this study are included in the manuscript and supporting files. Code is available on GitHub.
References

Evolution of dominance in metabolic pathwaysGenetics 168:1713–1735.https://doi.org/10.1534/genetics.104.028696

Effects of epistasis on phenotypic robustness in metabolic pathwaysMathematical Biosciences 184:27–51.https://doi.org/10.1016/S00255564(03)000579

Genetic interaction networks: toward an understanding of heritabilityAnnual Review of Genomics and Human Genetics 14:111–133.https://doi.org/10.1146/annurevgenom082509141730

Dynamic modeling of the central carbon metabolism of Escherichia coliBiotechnology and Bioengineering 79:53–73.https://doi.org/10.1002/bit.10288

Epistasis from functional dependence of fitness on underlying traitsProceedings of the Royal Society B: Biological Sciences 279:4156–4164.https://doi.org/10.1098/rspb.2012.1449

Bow ties, metabolism and diseaseTrends in Biotechnology 22:446–450.https://doi.org/10.1016/j.tibtech.2004.07.007

The synthetic genetic interaction spectrum of essential genesNature Genetics 37:1147–1152.https://doi.org/10.1038/ng1640

Empirical fitness landscapes and the predictability of evolutionNature Reviews Genetics 15:480–490.https://doi.org/10.1038/nrg3744

Missense meanderings in sequence space: a biophysical view of protein evolutionNature Reviews Genetics 6:678–687.https://doi.org/10.1038/nrg1672

Systematic mapping of genetic interaction networksAnnual Review of Genetics 43:601–625.https://doi.org/10.1146/annurev.genet.39.073003.114751

The causes and consequences of genetic interactions (epistasis)Annual Review of Genomics and Human Genetics 20:433–460.https://doi.org/10.1146/annurevgenom083118014857

The biomass objective functionCurrent Opinion in Microbiology 13:344–349.https://doi.org/10.1016/j.mib.2010.03.003

Systemic properties of metabolic networks lead to an epistasisbased model for heterosisTheoretical and Applied Genetics 120:463–473.https://doi.org/10.1007/s0012200912032

Evolution in the light of fitness landscape theoryTrends in Ecology & Evolution 34:69–82.https://doi.org/10.1016/j.tree.2018.10.009

Epistasis in a quantitative trait captured by a molecular model of transcription factor interactionsTheoretical Population Biology 77:1–5.https://doi.org/10.1016/j.tpb.2009.10.002

Epistasis and pleiotropy as natural properties of transcriptional regulationTheoretical Population Biology 49:58–89.https://doi.org/10.1006/tpbi.1996.0003

Modeling genetic architecture: a multilinear theory of gene interactionTheoretical Population Biology 59:61–86.https://doi.org/10.1006/tpbi.2000.1508