format: text

text:
# Abstract

## Abstract

Metabolic control analysis (Kacser & Burns (1973). Symp. Soc. Exp. Biol. 27, 65-104; Heinrich & Rapoport (1974). Eur. J. Biochem. 42, 89-95) was developed for the understanding of multi-enzyme networks. At the core of this approach is the flux summation theorem. This theorem implies that there is an invariant relationship between the control coefficients of enzymes in a pathway. One of the main conclusions that has been derived from the summation theorem is that phenotypic robustness to mutation ( e.g. dominance ) is an inherent property of metabolic systems and hence does not require an evolutionary explanation (Kacser & Burns (1981). Genetics. 97, 639-666; Porteous (1996). J. theor. Biol. 182, 223-232 ). Here we show that for mutations involving discrete changes ( of any magnitude ) in enzyme concentration the flux summation theorem does not hold. The scenarios we examine are two-enzyme pathways with a diffusion barrier, two enzyme pathways that allow for enzyme saturation and two enzyme pathways that have both saturable enzymes and a diffusion barrier. Our results are extendable to sequential pathways with any number of enzymes. The fact that the flux summation theorem cannot hold in sequential pathways casts serious doubts on the claim that robustness with respect to mutations is an inherent property of metabolic systems.

## 1. Introduction

How do changes in enzyme properties alter the physiological phenotype? Such knowledge is a key component in understanding the relation between genotype and phenotype. The two main approaches to this problem have developed into fields known as metabolic control analysis ( MCA ) \[Kacser & Burns, 1973, Heinrich & Rapoport, 1974, Kacser et al., 1995\] and biochemical systems theory ( BST ) \[Savageau, 1976, Savageau & Sorribas, 1989\]. Results from MCA have had extensive influence in biochemistry, genetics and evolution ( for an overview see the entire issue of J. theor. Biol., 182:3 , 1996 ). The cornerstone of the MCA approach is a theoretical result referred to as the flux summation theorem. The biological interpretation of the summation theorem is that there are systemic constraints inherent in metabolic pathways; these constraints limit the magnitude of the effects that changes in enzyme activity can have on flux through a pathway. The existence of such constraints has two implications. In the first place, the summation theorem implies that in general the control of flux in a pathway is shared between enzymes. Hence rate limiting enzymes are rare \[Kacser, 1995\]. Secondly, the theorem implies that on the average, mutations that change enzyme concentrations will have a small effect on flux. This implication has important consequences in evolutionary theory and genetics. It means that phenomena that fall under the rubric of phenotypic robustness to mutations ( such as selective neutrality, canalization and dominance ) are inherent properties of metabolic systems and are not a result of evolution. This MCA argument was originally made in reference to the case of dominance, stating that dominance is an inherent property of metabolic pathways and not a result of evolution \[Kacser & Burns, 1981, Keightley, 1996, Porteous, 1996\].

At the core of the MCA approach is a measure of the control exerted by each enzyme on flux through a pathway. The control coefficient $`C_i^J`$ measures how important each enzyme is in its ability to affect steady-state flux via changes in enzyme concentration. In its original formulation \[Kacser & Burns, 1973\] the control coefficient was defined as

$$C_i^J=\frac{\delta J}{J}/\frac{\delta E_i}{E_i},$$
(1)
where $`J`$ is the steady state flux of metabolites through the pathway (net rate of product formation at the end of a pathway), $`E_i`$ is the concentration of enzyme $`i`$, $`\delta E_i`$ is a discrete change in enzyme concentration and $`\delta J`$ is the resultant change in flux. Hence $`C_i^J`$ is basically a ratio of percentage change in a pathway’s steady state flux to the corresponding percentage change in enzyme concentration.

Theoretical results that contradict or object to different conclusions derived within MCA have been posed several times \[Cornish-Bowden, 1987, Giersch, 1988, Savageau & Sorribas, 1989, Savageau, 1992, Grossniklaus et al., 1996, Kholodenko et al., 1998\]. In each case the possible misrepresentation of non-linearities in MCA has appeared in a different guise. In the first place, Cornish-Bowden (1987) showed that in a sequential pathway if the maximal rate $`V_{max}`$ of consecutive enzymes are sequentially decreasing then dominance is not a necessary property of the pathway. Hence the possibility exists that dominance can evolve. This objection was rejected in MCA on grounds that the specific arrangement of kinetic values suggested is a special case that is very unlikely to occur by chance \[Kacser, 1987\]. A second set of objections was then put forth by Savageau and Sorribas (1989, 1992). They argued that in pathways exhibiting non-linear behavior that arises from properties such as enzyme-enzyme interactions, feedback loops or non-sequential pathway structure it can be shown that dominance is not inherent to the pathway. Based on disagreements on the mathematical models used, this second objection was rejected by MCA proponents \[Kacser, 1991, Kacser, 1995\]. Nevertheless more recent work indirectly corroborates this objection \[Grossniklaus et al., 1996, Kholodenko et al., 1998\].

The present paper deals with problems involving the biological in-applicability of MCA precepts. We show that for discrete changes ( of any magnitude ) in enzyme concentration, the flux summation theorem does not hold. This means that contrary to the MCA allegations, there are no a-priori constraints that would require the magnitude of mutational effects to be small. Furthermore these results re-open the possibility that phenotypic robustness to mutations ( e.g. dominance ) is not an inherent property of metabolic pathways and can be a result of evolution.

Our conclusions do not necessarily imply that the summation theorem is inherently false as an isolated mathematical proposition. Rather that it has been applied to the wrong context. This erroneous application arises due to two main disparities. In the first place, the argument for the summation theorem in the original Kacser and Burns paper (1973) is derived using a discrete formulation. However, in most subsequent work ( e.g. Kacser & Burns 1981 ) a continous formulation of the theorem is used. Hence the discrete proof for the summation theorem does not correspond to the continous form in which it is frequently used. In the second place, mutations necessarily involve discrete changes in enzyme activity. Hence, whether a continous version of the summation theorem holds or not ( which it would in cases where the flux function is homogenous, see Giersch 1988 ) does not say anything about whether a pathway will be robust or exhibit dominance with respect to discrete changes. In order for a pathway to be robust with respect to discrete changes, the discrete version of the summation theorem has to hold. In this paper we show that for dicrete changes of any size the discrete version does not hold, and hence that robustness is not an inherent property of a metabolic pathway.

The matter addressed here has been a subject of debate for thirty years. In order to be clear about the objects that we treat and the mathematical properites that pertain to them a rigorous approach to the problem was necessary. The main assumptions and implications we derived are summarized as follows:

#### INITIAL ASSUMPTION

The metabolic pathways studied are two-enzyme pathways. They have one input and one output and can reach a steady state flux.

#### CASE A : DIFFUSION LIMITED PATHWAYS

Maximal rate of input into the metabolic pathway is limited by a diffusion process.

The flux summation theorem does not hold.

#### CASE B : SEQUENTIAL PATHWAYS THAT ALLOW FOR ENZYME SATURATION

The two enzymes are near saturation.

The flux summation theorem does not hold.

#### CASE C : NUMERICAL ANALYSIS OF PATHWAYS WITH REVERSIBLE MICHAELIS-MENTEN KINETICS

Each enzyme in the pathway behaves according to reversible Michaelis-Menten enzyme kinetics. Isolated Michaelis-Menten enzymes do allow for enzyme saturation. Note that such kinetics are manifest in empirical studies of isolated enzymes in vitro.

When Michaelis-Menten enzymes are placed in the context of a multi-enzyme pathway, their kinetic mechanisms do not change. This is a pure assumption.

Maximal rate of input into the metabolic pathway is limited by a diffusion process. As in assumption 1, this is is an assumption that is empirically justified \[Tralau et al., 2000\].

In regimes where both enzymes are saturated the summation theorem does not hold.

In regimes where flux is near the diffusion limited rate, the summation theorem does not hold.

Mutations can cause heritable variation in the catalytic turnover rate $`k_{cat}`$ of an enzyme.

The control of a pathway can be switched from one enzyme to another via mutations affecting $`k_{cat}`$. Each enzyme can be the sole controlling enzyme of a pathway.

The main reason for the restricted applicability of the summation theorem is the fact that for its construction the non-linear interaction effects on flux were ignored. It is precisely the interaction effects between enzymes that lead to the manifestation of interaction effects ( epistasis ). Furthermore it is the interposition of the non-linear properties of enzymes that allow for the evolution of control.

## 2. Mathematical Introduction

### 2.1. The conception of control in MCA

Upon devising MCA, as Kacser and Burns put it, it was necessary to place the understanding of individual enzyme behavior into the context of a “metabolic society” composed of several enzymes \[Kacser & Burns, 1979\]. In accordance to this goal, one of the objectives of metabolic control analysis was to ascertain what role each enzyme could play in affecting steady state flux in a particular metabolic pathway. Specifically, how important each enzyme is in its ability to affect steady state flux via the regulation of the enzyme concentration. Hence the use of the control coefficient as $`C_i^J={\displaystyle \frac{\delta J}{J}}/{\displaystyle \frac{\delta E_i}{E_i}}`$, which is a non-dimensionalized sensitivity measure that shows the ratio of proportional change in flux to proportional change in enzyme concentration.

### 2.2. The flux summation theorem

One of the central tenets of metabolic control theory has been that the sum of the control coefficients in a pathway with $`n`$ enzymes equals one. That is,

$$\underset{i=1}{\overset{n}{}}C_i^J=1.$$
(2)
Equation (2) is commonly referred to as the flux summation theorem and was derived by Kacser and Burns using the following form of reasoning:

Given an unbranched chain of $`n`$ enzyme catalyzed reactions, if we were to increase the concentration $`E_i`$ of enzyme $`i`$ by the proportional ratio $`\delta E_i/E_i`$, the change in flux due to enzyme $`i`$ would be given by rearranging (1) to give:

$$\text{For all enzymes }i\text{ :}\left(\frac{\delta J}{J}\right)_i=\left(\frac{\delta E_i}{E_i}\right)C_i^J.$$
(3)
If $`\delta E_i/E_i=\alpha `$ , where $`\alpha `$ is a constant, one could consider the case where all enzyme concentrations are changed by the constant proportion $`\alpha `$. In such a case Kacser and Burns assume the following equation:

$$\frac{\delta J}{J}=\underset{i=1}{\overset{n}{}}\left(\frac{\delta J}{J}\right)_i.$$
(4)
Hence, using (4), (3) and $`\alpha =\delta E_i/E_i`$ ;

$$\frac{\delta J}{J}=\underset{i=1}{\overset{n}{}}\alpha C_i^J=\alpha \underset{i=1}{\overset{n}{}}C_i^J.$$
(5)
Separately, the assumption is made that if all enzyme concentrations are uniformly changed by $`\alpha `$ , then

$$\frac{\delta J}{J}=\alpha .$$
(6)
Substituting equation (6) into (5) we obtain equation (2):

$$\underset{i=1}{\overset{n}{}}C_i^J=1.$$

For a clearer view, the Kacser and Burns derivation of (2) can be summarized in the following conditionals:

$$\begin{array}{c}\left(\begin{array}{c}(\text{3});(\text{ Definition })\left[\left(\frac{\delta J}{J}\right)_i=\left(\frac{\delta E_i}{E_i}\right)C_i^J\right]\\ \\ (\text{4});(\text{ Assumption })\left[\frac{\delta J}{J}=_{i=1}^n\left(\frac{\delta J}{J}\right)_i\right]\\ \\ \text{( Restriction )}(i):[\delta E_i/E_i=\alpha ]\end{array}\right)\hfill \\ \\ \left(\begin{array}{c}\hfill (\text{5});\left[\frac{\delta J}{J}=\alpha _{i=1}^nC_i^J\right]\end{array}\right).\hfill \end{array}$$
(7)
$$\begin{array}{c}\left(\begin{array}{c}(\text{5});\left[\frac{\delta J}{J}=\alpha _{i=1}^nC_i^J\right]\\ \\ \text{( Restriction )}(i):[\delta E_i/E_i=\alpha ]\\ \\ (\text{6});(\text{ Assumption })[((i):[\delta E_i/E_i=\alpha ])[\frac{\delta J}{J}=\alpha ]]\end{array}\right)\hfill \\ \\ \left(\begin{array}{c}\hfill (\text{2});\left[_{i=1}^nC_i^J=1\right]\end{array}\right).\hfill \end{array}$$
(8)
Note that technically even if the assumptions (4) and (6) were true, (2) would be only proven for the restricted case where for every enzyme $`i`$ : $`[\delta E_i/E_i=\alpha ]`$.

## 3. Construction of two-enzyme pathway models

### 3.1. Two-enzyme pathways as dynamical systems

A general dynamical model of multi-enzyme systems can be formulated based on the classical model of single-substrate enzyme catalysis. Under such a scheme the existence of an intermediate enzyme-substrate complex is posited and a set of forward and backward reaction rates is attributed to each transformation. Consider a pathway in which an outside substrate $`u`$ diffuses into the initial reaction compartment,where the substrate is fed into a two reaction pathway that comprises two successive enzyme catalyzed reactions. An irreversible sink step is added to the end of the reaction sequence in order to remove the product.

$$\text{input(u)}\stackrel{\text{Diffusion}}{}s_1\stackrel{\text{Enzyme 1}}{}s_2\stackrel{\text{Enzyme 2}}{}s_3\stackrel{\text{Sink}}{}\text{output(o)}$$

The above scenario is represented by the following kinetic model shown schematically:

$$\left(\begin{array}{c}u\begin{array}{c}\stackrel{D}{}\\ \stackrel{D}{}\end{array}s_1\\ \\ e_1+s_1\begin{array}{c}\stackrel{k1}{}\\ \stackrel{k2}{}\end{array}(\underset{¯}{es}_1)\begin{array}{c}\stackrel{k3}{}\\ \stackrel{k4}{}\end{array}s_2+e_1\\ \\ e_2+s_2\begin{array}{c}\stackrel{k5}{}\\ \stackrel{k6}{}\end{array}(\underset{¯}{es}_2)\begin{array}{c}\stackrel{k7}{}\\ \stackrel{k8}{}\end{array}s_3+e_2\\ \\ s_3\stackrel{Q}{}o\end{array}\right)$$
(9)
where $`(\underset{¯}{es}_1)`$ and $`(\underset{¯}{es}_2)`$ are the enzyme-substrate complexes, while $`e1`$ and $`e2`$ are the first and second enzymes catalyzing the two respective reactions. The input and output correspond to $`u`$ and $`o`$ respectively. The kinetic model (9) is a class of model that can be represented as a system of autonomous differential equations with input $`u(t)`$ and initial conditions $`𝐬(0)`$ and $`𝐞(0)`$ such that;

$$\begin{array}{c}\dot{𝐬}(t)=f_1\text{[}𝐬(t),𝐞(t),u(t)\text{]}\\ \dot{𝐞}(t)=f_2\text{[}𝐬(t),𝐞(t)\text{]}\\ \dot{o}(t)=f_3\text{[}s_3(t)\text{]}\end{array}$$
(10)
where at time $`t`$, the vector $`𝐬(t)=s_1(t),s_2(t),s_3(t)`$ is the state vector of substrate concentrations, the vector $`𝐞(t)=e_1(t),\underset{¯}{es}_1(t),e_2(t),\underset{¯}{es}_2(t)`$ is the state vector of enzyme and enzyme-substrate complex concentrations, the variable $`o(t)`$ is the concentration associated with the output variable and $`\dot{o}(t)`$ is the output flux. The notation $`\dot{x}`$ associated with any variable $`x`$ denotes the time derivative $`dx/dt`$. If the system reaches a steady state flux the output flux $`\dot{o}(t)`$ is equivalent to the steady state flux $`J`$ used in the MCA literature.

### 3.2. Characterization of two-enzyme pathways into distinct classes

We will make the rest of our argument using a set theroetic framework. For any given set, if we can define the general properties of its members, then any object which is deemed to be a member of that set will also have those properties. A system of equations is a logical proposition. As such we can cast our dynamical system model in a set theoretic guise. A pathway $`p`$ is an object. Each object has a series of properties $`K`$ associated with it, such that $`Kp=\text{True}`$. The properties $`K`$ delineate the context and system of equations that can describe a pathway ( For alternative representations see Fontana & Buss 1996 ). Each pathway $`p`$ can then be treated as an element of a set with distinct properites.

#### Definition 1 : $`𝐏_{\mathrm{𝟏},\mathrm{𝟏}}`$ $``$ Set of two enzyme pathways with one input and one output.

There are no strict stipulations on the stoichiometry nor on the specifics of enzyme kinetic mechanisms for membership in $`𝐏_{\mathrm{𝟏},\mathrm{𝟏}}`$ ( see Appendix A ).

#### Definition 2: $`𝐏_{𝐮,𝐉}`$ Set of two-enzyme pathways with constant input and steady-state flux

The set $`𝐏_{𝐮,𝐉}`$ is defined as the subset $`𝐏_{𝐮,𝐉}𝐏_{\mathrm{𝟏},\mathrm{𝟏}}`$ such that $`p𝐏_{𝐮,𝐉}`$ if and only if the following four conditions hold:

1. $`p𝐏_{\mathrm{𝟏},\mathrm{𝟏}}`$
2. Irrespective of the enzyme kinetic mechanisms, $`p`$ has the stoichiometry of a two-enzyme sequential pathway as exhibited in model (9).
3. When driven by a constant input $`u`$, $`p`$ exhibits a unique steady state output flux $`J^{}`$ that can be associated with each pair of total enzyme concentrations $`E_1,E_2`$.
4. For any error margin $`\epsilon >0`$, the steady state output flux $`J^{}`$ can be approximated by using the total enzyme concentrations $`E_1,E_2`$ and constant input $`u`$ as arguments for a computable function $`g_\epsilon `$ such that $`J=g_\epsilon [u,E_1,E_2]`$ and $`|JJ^{}|<\epsilon `$ .

Items 2 to 4 are elaborated more precisely in Definitions 2.2 to 2.4 in Appendix A. Note that the MCA approach holds informal versions of these propositions as premises, without which talking about the steady-state flux $`J`$ and the summation theorem would not be possible.

#### Definition 3: $`𝐏_{\mathrm{𝐮𝐃},𝐉}`$ Set of diffusion limited two-enzyme pathways with steady-state flux

The set $`𝐏_{\mathrm{𝐮𝐃},𝐉}`$ is defined as the subset $`𝐏_{\mathrm{𝐮𝐃},𝐉}𝐏_{𝐮,𝐉}`$ such that $`p𝐏_{\mathrm{𝐮𝐃},𝐉}`$ if and only if the following two conditions hold:

1. $`p𝐏_{𝐮,𝐉}`$
2. The net input flux from $`u(t)`$ into the system is governed by a diffusion-limited process such that $`Netstartflux=D\left(u(t)s_1(t)\right)`$ where $`D`$ is a diffusion constant.

Item 2 is more precisely defined in definition 3.2 in Appendix A.

#### Definition 4: $`𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ Set of two-enzyme pathways with steady-state flux and kinetics that allows for enzyme saturation.

The set $`𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ is defined as the subset $`𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}𝐏_{𝐮,𝐉}`$ such that $`p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ if and only if the following two conditions hold:

1. $`p𝐏_{𝐮,𝐉}`$
2. For every enzyme $`i`$, there exists a constant $`k_{cat(i)}0`$ such that for all values of $`E_i0`$ the steady state flux is consistent with the relation $`JE_ik_{cat(i)}`$.

Item 2 is more precisely defined in definition 4.2 in Appendix A.

#### Definition 5: $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ Set of diffusion limited two-enzyme pathways with steady-state flux and kinetics that allows for enzyme saturation.

The set $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ is defined as the intersection

$$𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}=(𝐏_{\mathrm{𝐮𝐃},𝐉}𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉})𝐏_{𝐮,𝐉}.$$

As an intuitive guide, the order of the generality of the representations discussed so far, in order of higher generality are: $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}𝐏_{\mathrm{𝐮𝐃},𝐉}𝐏_{𝐮,𝐉}𝐏_{\mathrm{𝟏},\mathrm{𝟏}}`$ and $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}𝐏_{𝐮,𝐉}𝐏_{\mathrm{𝟏},\mathrm{𝟏}}`$ (See Figure 1).

## 4. Problems with the summation theroem

The definitions from the previous section have a series of logical implications. Here we present these implications and explain their significance. The proofs for these propositions are given in Appendix B.

### 4.1. Pathways with a constant input and a steady-state flux $`(𝐏_{𝐮,𝐉})`$

#### Proposition 11

For any pathway of class $`𝐏_{𝐮,𝐉}`$, the flux summation theorem implies that for any given ratio of enzymes $`0(E_1/E_2)\mathrm{}`$ , there exists a $`C_\theta `$ such that flux is a linear function of the magnitude of the vector $`E_1,E_2`$.

$$\begin{array}{c}(p𝐏_{𝐮,𝐉}):(\theta |\mathrm{tan}\theta =E_1/E_2):(C_\theta ):\hfill \\ \left[C_1^J+C_2^J=1J=E_1,E_2C_\theta \right].\hfill \end{array}$$
(11)
The meaning of proposition 11 is very simple. If a pathway were to be of class $`𝐏_{𝐮,𝐉}`$ and were to obey the flux summation theorem the following would hold: the flux would behave in a linear fashion with respect to tandem changes in which both enzyme concentrations are changed by the same percentage. It is a short few steps (Appendix B) to show that

$$(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}):[\frac{\delta J}{J}=\alpha J=E_1,E_2C_\theta ].$$
(12)
Proposition (12) is equivalent to a statement of homogeneity ( See Giersch 1988 ), such that if $`J=g[E_1,E_2]`$ then

$$g[(1+\alpha )E_1,(1+\alpha )E_2]=(1+\alpha )g[E_1,E_2].$$
(13)
For example consider two enzymes whose genes are located on the same operon. In addition consider the case in which such a cell is in a state in which the concentrations of the two enzymes are initially negligible. When the operon is turned on, each enzyme will have a particular rate of expression. If the ratio between those rates of expression remains constant, then the flux summation theorem would imply that the flux would increase linearly without reaching a saturation plateau. This is what $`J=E_1,E_2C_\theta `$ implies, where $`E_1,E_2=\sqrt{E_1^2+E_2^2}`$. Hence, in a pathway of class $`𝐏_{𝐮,𝐉}`$, in which we have made no assertion about the particulars of the enzyme-kinetic equations or parameters, the flux summation theorem assumes a particular type of dynamic behavior. A flux that increases linearly without saturation is a very special case and such a condition can only hold in a restricted subset of $`𝐏_{𝐮,𝐉}`$. In fact, as proposition 14 will attest, even for a simple pathway with the sole kinetic imposition of a diffusion barrier, the flux summation theorem fails.

### 4.2. Pathways with a diffusion barrier $`(𝐏_{\mathrm{𝐮𝐃},𝐉})`$

#### Proposition 14

For any pathway of class $`𝐏_{\mathrm{𝐮𝐃},𝐉}`$, the flux summation theorem is false.

$$(p𝐏_{\mathrm{𝐮𝐃},𝐉}):[(\delta E_1,\delta E_2):(E_1,E_2_0^1):[C_1^J+C_2^J1]]$$
(14)
Proposition 14 basically states the general result that any pathway which has a diffusion barrier at a point before its first step will not obey the summation theorem. This includes almost any metabolic pathway that takes place in the confines of a cell. In fact it can be shown that any pathway that exhibits a plateauing effect on its flux surface will not obey the summation theorem. The proof of this theorem is intuitively simple (Appendix B). It is basically dependent on the conditional relation

$$C_1^J+C_2^J=1J=E_1,E_2C_\theta $$
(15)
from (11). For any pathway class that is a subset of $`𝐏_{𝐮,𝐉}`$, if the pathways do not obey the right hand side of 15 then they do not obey the left hand side either. For $`𝐏_{\mathrm{𝐮𝐃},𝐉}`$ one can show that $`J=E_1,E_2C_\theta `$ does not hold when $`C_\theta >0`$. This is because the diffusion barrier will impede flux from increasing indefinitley. This implies that $`_{i=1}^nC_i^J=1`$ cannot be true for all enzyme concentrations. This falsehood applies to all diffusion-limited flux surfaces except for a disfunctional cell with flux surface $`J=0`$. Proposition 14 is due to the fact that control is not only dependent on the enzymatic steps in a pathway, but also the diffusion steps.

### 4.3. Pathways that allow for enzyme saturation $`(𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉})`$

#### Proposition 16

For any pathway of class $`𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$, if both enzymes approach saturation and changes in enzyme concentration are positive, then the sum of control coefficients approaches zero.

$$(p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}):[(\delta E_1,\delta E_2_{>0}^1):\underset{Sat_1,Sat_21}{lim}(C_1^J+C_2^J)=0].$$
(16)
Whenever two enzymes are near saturation, a situation is created whereupon the flux effects of increasing the concentration for one enzyme is restrained by the other enzyme. Hence whenever two enzymes approach saturation simultaneously, the sum of their control coefficients approaches zero and the flux summation theorem does not hold ( see appendix B ). In fact it is not too difficult to extend the proof such that for a sequential pathway of any length, the sum of all control coefficients will approach zero whenever any two enzymes approach saturation.

#### Proposition 17

For any pathway of class $`𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ in which a decrease in enzyme concentration increases saturation, if both enzymes approach saturation and changes in enzyme concentration are negative, then the sum of control coefficients approaches two.
For all intervals $`a_i<E_i<b_i`$:
If $`p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$ and $`i`$ : ( $`\frac{Sat_i}{E_i}0`$ and $`a_i0`$ and $`\delta E_i<0`$ ) then:

$$\underset{Sat_1,Sat_21}{lim}\left(C_1^J+C_2^J\right)=2.$$
(17)
In such situations the pathway is no longer robust to decreases in enzyme concentration. Hence dominance is not an inherent property of these pathways. It is simple to extend this result to sequential pathways of any length, in which the sum of control coefficients will approach the the number of enzymes which are approaching saturation.

From propositions 14, 16 and 16 it also follows that the summation theorem is also false in pathways of class $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$.

### 4.5. Enzyme interactions in $`𝐏_{𝐮,𝐉}`$

The motivation for analyzing multi-enzyme pathways is that they behave differently from single-enzymes. The natural question is to ask what are the direct or indirect interactions between enzymes that lead to the difference between single-enzyme and multi-enzyme pathway behavior. Interaction does not necessarily have to signify direct physical interaction between two enzymes. It can occur via the effect of one enzyme on the concentration of pathway intermediates that in turn physically interact with another enzyme. In equation (6) which was one of the assumptions that Kacser and Burns made in deriving the summation theorem, the proposition $`\frac{\delta J}{J}=_{i=1}^n\left(\frac{\delta J}{J}\right)_i`$ amounts to an assumption that enzymes do not exhibit interactions. Here interaction is specifically referring to the mutual effect that the two enzymes would have on each-other’s ability to affect pathway flux:

#### Proposition 18

For any pathway of class $`𝐏_{𝐮,𝐉}`$, the proposition $`\left[\frac{\delta J}{J}=_{i=1}^n\left(\frac{\delta J}{J}\right)_i\right]`$ (4) and the restriction $`\left(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}\right)`$ imply that the change in flux due to change of total enzyme concentration at one locus is independent of total enzyme concentration at the second locus in the pathway.

$$\begin{array}{c}(p𝐏_{𝐮,𝐉}):(E_1,E_2_0^1):\hfill \\ [[(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}):\frac{\delta J}{J}=_{i=1}^2\left(\frac{\delta J}{J}\right)_i]\frac{^2J}{E_1E_2}=0].\hfill \end{array}$$
(18)
The assumption of independence of effects is something that surreptitiously goes against the motivation for developing MCA. Nonetheless the assumption of no enzyme interactions is built into the summation theorem from the beginning. In fact:

#### Proposition 19

For any pathway of class $`𝐏_{𝐮,𝐉}`$, the summation theorem implies that the change in flux due to change of total enzyme concentration at one locus is independent of total enzyme concentration at the second locus in the pathway.

$$\begin{array}{c}(p𝐏_{𝐮,𝐉}):(E_1,E_2_0^1):\hfill \\ [[(\delta E_1,\delta E_2):C_1^J+C_2^J=1]\frac{^2J}{E_1E_2}=0].\hfill \end{array}$$
(19)
This means that the summation theorem implies the inexistence of epistasis. It immediately follows that

$$\begin{array}{c}(p𝐏_{\mathrm{𝐮𝐃},𝐉}):\left(E_1,E_2_0^1|\frac{^2J}{E_1E_2}0\right):\hfill \\ \left[(\delta E_1,\delta E_2):\left[C_1^J+C_2^J1\right]\right].\hfill \end{array}$$
(20)
Hence in any regime where epistasis is present the summation theorem is violated. Furthermore:

#### Proposition 21

For any pathway of class $`𝐏_{𝐮,𝐉}`$, the summation theorem holds if and only if flux is a linear function of enzyme concentrations.

$$\begin{array}{c}(p𝐏_{𝐮,𝐉}):(E_1,E_2_0^1):\hfill \\ [[(\delta E_1,\delta E_2):C_1^J+C_2^J=1]J=E_1c_1+E_2c_2].\hfill \end{array}$$
(21)
In short, the summation theorem imposes conditions on enzyme interactions that are neither explicitly argued for nor conceptually justified.

## 5. Numerical analysis of a pathway of class $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$

Here we use numerical simulations for a particular case in order to verify the analytical results from the previous section. Consider a case where the enzymes exhibit reversible Michaelis-Menten kinetics. A diffusion barrier is included at the beginning of the pathway. During the interval of analysis $`0tt_f`$, the input remains constant such that $`\dot{\overline{u}}(t)=0`$ (e.g. approximation of large nutrient pool and/or steady state conditions of a chemostat). In such a case the full form of the system of differential equations (10) combined with the kinetic model (9) are:

$$\left(\begin{array}{c}\frac{d\left[s_1\right]}{dt}=D[s_1]k1[e_1][s_1]+D[u]+k2[\underset{¯}{es}_1]\hfill \\ \\ \frac{d\left[\underset{¯}{es}_1\right]}{dt}=k1e_1[s_1]k2[\underset{¯}{es}_1]k3[\underset{¯}{es}_1]+k4[e_1][s_2]\hfill \\ \\ \frac{d\left[s_2\right]}{dt}=k3[\underset{¯}{es}_1]k4[e_1][s_2]k5[e_2][s_2]+k6[\underset{¯}{es}_2]\hfill \\ \\ \frac{d\left[\underset{¯}{es}_2\right]}{dt}=k5[e_2][s_2]k6[\underset{¯}{es}_2]k7[\underset{¯}{es}_2]+k8[e_2][s_3]\hfill \\ \\ \frac{d\left[s_3\right]}{dt}=k7[\underset{¯}{es}_2]Q[s_3]k8[e_2][s_3]\hfill \\ \\ \frac{d\left[e_1\right]}{dt}=k1[e_1][s_1]+k2[\underset{¯}{es}_1]+k3[\underset{¯}{es}_1]k4[e_1][s_2]\hfill \\ \\ \frac{d\left[e_2\right]}{dt}=k5[e_2][s_2]+k6[\underset{¯}{es}_2]+k7[\underset{¯}{es}_2]k8[e_2][s_3]\hfill \\ \\ \frac{d\left[o\right]}{dt}=Q[s_3]\hfill \end{array}\right)$$
(22)
where the postfix $`(t)`$ for each state variable is omitted.

The system of differential equations (22) will serve as the basis for the examination of a particular case. For the present analysis we are interested in the steady states of the dynamical system. At steady state, as defined in appendix A, the derivatives $`\dot{𝐬}(t)\text{ and }\dot{𝐞}(t)`$ are equal to zero. Consider a case where the equations in (22) serve as the basis for the functions $`f_1`$, $`f_2`$ and $`f_3`$ with the stoichiometric constraint that

$$\begin{array}{c}[E_1]=[e_1]+[\underset{¯}{es}_1]\\ [E_2]=[e_2]+[\underset{¯}{es}_2].\end{array}$$
(23)
In such a case, the steady state solution can be obtained in an implicit form ( See appendix C ). The implicit equations can then be solved numerically by Newton-Raphson based algorithms that converge to the flux surface shown in Figure 2. On a cursory level the surface shows the plateau effect familiar in most models of two-enzyme pathways. The plateau arises due to the limiting effects of the diffusion barrier. We will refer to the pathway based on equations (22) as pathway $`p_{\text{22}}`$. It is not difficult to discern that pathway $`p_{\text{22}}`$ is a member of the set $`𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$. It satisfies the conditions required in definitions 1 to 5.

### 5.1. Test of the summation theorem

Given the flux surface on $`p_{\text{22}}𝐏_{\mathrm{𝐮𝐃},𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$, the flux summation theorem can be tested directly. Figures 3 and 4 show surfaces for for discrete versions of the sum $`_{i=1}^nC_i^J`$ for $`p_{\text{22}}`$. It is evident from these examples that in the pathway $`p_{\text{22}}`$ there is no evidence of an invariant relation in which the control coefficients sum to the predetermined value of one. Figures 3 and 4 serve as illustrations of propositions 14, 16 and 17. The summation theorem fails in the regimes where flux is near its diffusion-limited maximal rate or in areas where both enzymes are near saturation. The regions where both enzymes approach saturation are shown in Figure 5.

### 5.2. Epistatic interactions between enzymes

The question of the extent to which enzymes have interactive effects on flux can also be addressed directly using the pathway $`p_{\text{22}}`$ as an example. In Figure 6 the mixed partial difference $`\delta ^2J/\delta E_1\delta E_2`$ is used as a measure of epistatic interactions between enzymes (proposition 18). Note that the region where the epistatic interactions increase coincides with a region where the summation theorem fails and both enzymes are approaching saturation.

### 5.3. Rate limiting enzymes and the switching of control between enzymes

One of the claims from MCA has been that there is no such thing as a rate limiting enzyme in a pathway. The associated claim is that control is mostly shared between enzymes in a pathway \[Kacser, 1995\]. The results of our numerical simulations of $`p_{\text{22}}`$ do not support this claim. Figure 7 shows the numerical results for two pathways with different $`k_{cat}`$ configurations. In the range of concentration regimes $`E_1,E_2`$ examined, most regimes manifest a sole enzyme as the controlling enzyme for a given pathway. In regimes where flux is near the diffusion limited maximal flux, no enzyme has control. The only region where the enzymes share control is where the summation theorem fails. In fact in the regions that the summation theorem holds, the sum $`_{i=1}^2C_i^J=1`$ is basically the value $`C_z^J`$ for the sole controlling enzyme $`z`$ at that regime. Furthermore the identity of the controlling enzyme at a given regime is not invariant. The identity of the controlling enzyme can change by changing the $`k_{cat}`$ values of the enzymes involved (e.g. by mutation ). In Figure 7 the switch in control is clearly apparent by comparing the values of $`C_2^J`$ between the two $`k_{cat}`$ configurations.

## 6. Discussion and conclusion

The main intent of the present paper has been to present an analysis and critique of some of the main biological conclusions derived from MCA. Our results strongly support initial objections that were made with regards to the MCA perspective on the biochemical nature of dominance ( Cornish-Bowden 1987, Savageau 1992. See also Mayo & Burger 1997 ). At the core of the problem of multienzyme systems is interaction effects between enzymes. In order to address these interactions, the non-linear properties of enzyme catalysis cannot be ignored. These non-linear interactions play a major role in the phenotypic characteristics of metabolic physiology. Our analysis shows that for discrete changes of any magnitude, the flux summation theorem does not hold. This means that there are no a-priori constraints that would require the magnitude of mutational effects to be in a low range. This leaves open the possibility for two phenomena which had been previously rejected in MCA. In the first place, rate limiting steps or controlling enzymes do not have to be rare, and their identity can change by evolution. Secondly, phenotypic robustness with respect to mutations is not an automatic property of metabolic pathways; the degree of robustness can be modified through the evolution of kinetic properties. These are implications that we will examine in a later work.

## Appendices

## A. Definitions

In this section some of the conditions proposed in definitions 1,2, 3 and 4 are formulated in a more precise language as definitions 2.2, 2.3, 2.4, 3.2 and 4.2. In these definitions quantifiers are encapsulated in parenthesis and propositions in square brackets. The symbol “$``$” denotes the boolean expression “and”. The symbol “$`_0^1`$” denotes the real numbers greater than or equal to zero. Ordered tuples and vectors are encapsulated by the symbols “$``$”. The property $`K`$ associated with an object $`p`$ is denoted as $`(\text{“name of property”})p`$.

### Definition 1

$$\begin{array}{c}\\ (p𝐏_{\mathrm{𝟏},\mathrm{𝟏}}):\begin{array}{c}p\left[\begin{array}{c}\dot{𝐬}(t)=f_1[𝐬(t),𝐞(t),u(t)]\hfill \\ \dot{𝐞}(t)=f_2[𝐬(t),𝐞(t)]\hfill \\ \dot{o}(t)=f_3[s_3(t)]\hfill \\ \left(u(t),o(t)_0^1\right)\left(𝐞(0)_0^4\right)\left(𝐬(0)_0^3\right)\hfill \end{array}\right]\hfill \end{array}\hfill \end{array}$$
(24)
### Definition 2.2

The proposition $`(\mathrm{`}\mathrm{`}\text{Stoichiometry: two enz. seq.}\mathrm{"})p`$ imposes both organizational and kinetic constraints on the pathway $`p𝐏_{𝐮,𝐉}`$ by designating the stoichiometric properties of the pathway. Such a proposition is necessary for an abstract representation of physical mass balance.

$$\begin{array}{c}(\mathrm{`}\mathrm{`}\text{Stoichiometry: two enz. seq.}\mathrm{"})p\hfill \\ \begin{array}{c}((t,u(t),o(t))_0^1):(𝐞(0)_0^4):(𝐬(0)_0^3):\hfill \\ \left(\begin{array}{c}\left[\begin{array}{c}\begin{array}{c}[\dot{E}_1(t)=\dot{E}_2(t)=0]\hfill \\ [E_1=E_1(t)=e_1(t)+\underset{¯}{es}_1(t)]\hfill \\ [E_2=E_2(t)=e_2(t)+\underset{¯}{es}_2(t)]\hfill \\ [E_10E_20]\hfill \end{array}\hfill \end{array}\right]\hfill \\ \\ \left(\begin{array}{c}\left[\begin{array}{c}\xi (t)=s_1(t)+2\underset{¯}{es}_1(t)+s_2(t)+2\underset{¯}{es}_2(t)\hfill \\ +s_3(t)+e_1(t)+e_2(t)\hfill \end{array}\right]\hfill \\ \\ \left(\begin{array}{c}(I_{in}(t),I_{rev}(t)_0^1):\hfill \\ [\dot{\xi }(t)=I_{in}(t)I_{rev}(t)\dot{o}(t)]\hfill \end{array}\right)\hfill \end{array}\right)\hfill \end{array}\right)\hfill \end{array}\hfill \end{array}$$
(25)
where the pair $`E_1,E_2_0^2`$ are the total enzyme concentrations for any pathway $`p`$ such that $`E_1,E_2`$ does not change for each particular vector of starting conditions and that for a given set of initial conditions, the concentrations $`E_1\text{ and }E_2`$ are to be determined at any time $`t`$ by summing the free-enzyme and enzyme-substrate complex concentrations for each enzyme respectively. The variable $`\xi (t)`$ needs some explanation. $`\xi (t)`$ determines the “stoichiometric content” of the pathway in the sense that the coefficients $`c_i`$ in the equation

$$\xi (t)=c_1s_1(t)+c_2\underset{¯}{es}_1(t)+c_3s_2(t)+c_4\underset{¯}{es}_2(t)+c_5s_3(t)+c_6e_1(t)+c_7e_2(t)$$

are determined by whether in reference to the input molecule, the transformations leading to the macro-molecules associated with each coefficient are fusions, cleavages or one-one transformations. The values of the $`c_i`$ in (25) correspond to the kinetic model (9). The equation for $`\dot{\xi }(t)`$ also deserves attention. Each metabolic pathway can be characterized such that $`\dot{\xi }(t)`$ is determined by the additive combination of two values at time $`t`$:

$$\dot{\xi }(t)=\text{Netstartflux(t) - Netendflux(t)}$$

which in turn can be subdivided as

$$\begin{array}{c}\text{Netstartflux(t) = }\mathrm{I}_{\mathrm{i}\mathrm{n}}(\mathrm{t})\mathrm{I}_{\mathrm{r}\mathrm{e}\mathrm{v}}(\mathrm{t})\text{ }\text{ and}\hfill \\ \text{Netendflux(t) =}O_{out}(t)O_{rev}(t)\hfill \end{array}$$

where $`I_{in}(t)`$ is the flux entering the system from the input variable $`u(t)`$, $`I_{rev}(t)`$ is the flux exiting the system and going back into $`u(t)`$, $`O_{out}(t)`$ is the flux exiting the system into the output variable $`o(t)`$ and $`O_{rev}(t)`$ is the flux entering back into the system from $`o(t)`$. For a kinetic model with an irreversible sink step, $`O_{rev}(t)=0`$. For a system where the output variable $`o(t)`$ is only driven by a single pathway $`p`$, $`O_{out}(t)=\dot{o}(t)`$. Finally for $`𝐏_{𝐮,𝐉}`$ the only stipulations made for the two components of Netstartflux(t) is that they must both be greater or equal to zero such that $`I_{in}(t),I_{rev}(t)_0^1`$ and $`\text{Netstartflux(t)}=I_{in}(t)I_{rev}(t)`$.

### Definition 2.3

For a pathway $`p𝐏_{𝐮,𝐉}`$ the logical proposition
$`(\mathrm{`}\mathrm{`}\text{Const. u }\text{ Unique }u,E_1,E_2,J^{}\mathrm{"})p`$ states that given an input $`u(t)`$, if $`u(t)`$ remains constant then the pathway will reach a unique steady state flux $`J^{}_0^1`$ that is dependent on total enzyme concentrations $`E_1\text{ and }E_2`$ such that:

$$\begin{array}{c}(\mathrm{`}\mathrm{`}\text{Const. u }\text{ Unique }u,E_1,E_2,J^{}\mathrm{"})p\hfill \\ (u(0)_0^1):(E_1(0),E_2(0)_0^1):\hfill \\ (J^{}_0^1):(o(t)_0^1):(s_1(0),s_2(0),s_3(0)_0^1):\hfill \\ \left(\begin{array}{c}\left(\begin{array}{c}(t_0^1):\hfill \\ [\dot{u}(t)=0]\hfill \end{array}\right)\left[\begin{array}{c}\underset{t\mathrm{}}{lim}\dot{𝐬}(t)=0,0,0\hfill \\ \underset{t\mathrm{}}{lim}\dot{𝐞}(t)=0,0,0,0\hfill \\ \underset{t\mathrm{}}{lim}\ddot{o}(t)=0\hfill \\ \underset{t\mathrm{}}{lim}\dot{o}(t)=J^{}\hfill \end{array}\right]\hfill \end{array}\right)\hfill \end{array}$$
(26)
Note that at steady state $`\dot{\xi }(t)=0`$.

### Definition 2.4

The proposition $`(\mathrm{`}\mathrm{`}\text{computable: }E_1,E_2,uJ\mathrm{"})p`$ states that an approximation $`J`$ of $`J^{}`$ is computable as a function of the initial input and total enzyme concentrations:

$$\begin{array}{c}(\mathrm{`}\mathrm{`}\text{computable: }E_1,E_2,uJ\mathrm{"})p\hfill \\ \left(\begin{array}{c}(\epsilon _{>0}^1):(g_\epsilon :_0^3^1):\hfill \\ ((u(0),o(t))_0^1):(𝐞(0)_0^4):(𝐬(0)_0^3):\hfill \\ \left(\begin{array}{c}[J=g_\epsilon [u(0),E_1,E_2]|JJ^{}|\epsilon ]\hfill \end{array}\right)\hfill \end{array}\right)\hfill \end{array}$$
(27)
The above definition means that for any error margin $`\epsilon 0`$, there exists at-least one effectively computable function $`g_\epsilon `$ that maps to the approximate value $`J`$ by using the total enzyme concentrations and initial input as arguments. Note that $`g_\epsilon `$ as an effectively computable function refers to both numerical and analytic functions that can produce an approximation $`J`$. The assumption of computability is associated with the assumption that there is a rule based mechanistic relation between $`E_1`$, $`E_2`$ and $`u`$. There is a strong assumption inherent in the definitions of $`g_\epsilon `$ in $`𝐏_{𝐮,𝐉}`$: that for any pathway $`p𝐏_{𝐮,𝐉}`$ a steady state flux exists and can be determined using only $`u(0)`$, $`E_1`$ and $`E_2`$ as arguments. That is, except for the constraints imposed by $`u(0),E_1`$ and $`E_2`$ the steady state is independent of the specifics of starting conditions. Note that for classes of pathways that exhibit multiple equilibria the trajectories become important and all the starting conditions have to be included as arguments in a newly defined computable function other than $`g_\epsilon `$, which would be embedded in a system other than $`𝐏_{𝐮,𝐉}`$.

### Definition 3.2

$$\begin{array}{c}(\mathrm{`}\mathrm{`}\text{Diffusion Limited }\text{I}\mathrm{"})p\hfill \\ (D_0^1):\hfill \\ ((t,u(t),o(t))_0^1):(𝐞(\mathrm{𝟎})_\mathrm{𝟎}^\mathrm{𝟒}):(𝐬(\mathrm{𝟎})_\mathrm{𝟎}^\mathrm{𝟑}):\hfill \\ \left(\begin{array}{c}[I_{in}(t)=Du(t)]\hfill \\ [I_{rev}(t)=Ds_1(t)]\hfill \end{array}\right)\hfill \end{array}$$
(28)
where $`D`$ is the diffusion constant.

### Definition 4.2

$$\begin{array}{c}(\mathrm{`}\mathrm{`}\text{Saturable Enzymes}\mathrm{"})p\hfill \\ ((k_{cat(1)},k_{cat(2)})_0^1):((E_1,E_2)_0^1):\hfill \\ ((t,u(t),o(t))_0^1):(𝐞(\mathrm{𝟎})_\mathrm{𝟎}^\mathrm{𝟒}):(𝐬(\mathrm{𝟎})_\mathrm{𝟎}^\mathrm{𝟑}):\hfill \\ \left(\begin{array}{c}[JE_1k_{cat(1)}]\hfill \\ [JE_2k_{cat(2)}]\hfill \end{array}\right)\hfill \end{array}$$
(29)
## B. Proofs

### B.1 Proof of proposition 11

We start with the summation theorem $`_{i=1}^nC_i^J=1`$. Using the definition of $`C_i^J`$ from (1) we have

$$(E_1,E_2_{>0}^1):\frac{\delta J}{\delta E_1}\frac{E_1}{J}+\frac{\delta J}{\delta E_2}\frac{E_2}{J}=1.$$
(30)
The continuous version of (30) is $`\underset{\delta E_1,\delta E_20}{lim}\left(\frac{\delta J}{\delta E_1}E_1+\frac{\delta J}{\delta E_2}E_2\right)=J`$ which is equivalent to

$$\frac{J}{E_1}E_1+\frac{J}{E_2}E_2=J.$$
(31)
Equation (31) is a first order quasi-linear partial differential equation with general solution

$$J=E_2f\left[\frac{E_1}{E_2}\right]=E_1f^{}\left[\frac{E_2}{E_1}\right],$$
(32)
where $`f,f^{}:^1^1`$ (see for example Rade & Westergren 1995, p230). For each fixed ratio of $`E_1/E_2`$, the vector $`E_1,E_2`$ forms a fixed angle $`\theta `$ with the coordinate axes $`E_2`$ and $`E_1`$ such that $`E_1/E_2=\mathrm{tan}\theta `$. Hence we can formulate the following proposition:

$$\begin{array}{c}(g_{st}:_0^2^1):(E_1,E_2_{>0}^1):(f_{st}:_0^1^1):\\ J=g_{st}[E_1,E_2]=E_2f_{st}[\mathrm{tan}\theta ]=h\mathrm{cos}\theta f_{st}[\mathrm{tan}\theta ]\end{array}$$
(33)
where $`h`$ is the norm of the vector $`E_1,E_2`$ such that $`h=E_1,E_2=\sqrt{E_1^2+E_2^2}`$. Given that for every $`\theta `$ there exists a $`C_\theta `$ such that $`C_\theta =\mathrm{cos}\theta f_{st}[\mathrm{tan}\theta ]`$, we have

$$(\theta 0\theta \pi /2):(C_\theta ):[J=hC_\theta ].$$
(34)
$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{11}`$

### B.2 Proof of Proposition 12

The proof of the proposition $`\left[\frac{\delta J}{J}=\alpha J=E_1,E_2C_\theta \right]`$ is uncomplicated. We start with proving the forward conditional $`\left[\frac{\delta J}{J}=\alpha J=E_1,E_2C_\theta \right]`$. By definition we are given that

$$\alpha =\delta E_1/E_1=\delta E_2/E_2.$$
(35)
Consider the vector $`E_1,E_2`$ and its norm $`h=E_1,E_2=\sqrt{E_1^2+E_2^2}`$. Given simultaneous increases in enzyme concentration $`\delta E_1`$ and $`\delta E_2`$ that obey (35) we have $`h+\delta h=\sqrt{(1+\alpha )^2E_1^2+(1+\alpha )^2E_2^2}`$. Hence $`\delta h/h=\alpha `$ and thereby

$$\delta J/J=\delta h/h.$$
(36)
In its continuous form equation (36) can be set up as a separable differential equation such that

$$(E_1/E_2):\left[\frac{dJ}{J}=\frac{dh}{h}\right].$$
(37)
Integrating both sides of (37) we have the solution

$$(E_1/E_2):(C_\theta ):\left[J=hC_\theta \right].$$
(38)
The backward conditional $`\left[\frac{\delta J}{J}=\alpha J=E_1,E_2C_\theta \right]`$ can be proven by differentiating (38) with respect to $`h`$.

$`\mathrm{}\mathrm{𝐏𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{12}`$

### B.3 Proof of proposition 14

Consider a case where we posit the existence of a pathway $`p_1𝐏_{𝐮,𝐉}`$ which obeys the flux summation theorem . In such a case, using (34) we have;

$$\begin{array}{c}(\theta 0\theta \pi /2):\hfill \\ [[h\mathrm{}][E_1\mathrm{}E_2\mathrm{}]].\hfill \end{array}$$
(39)
Furthermore for the function $`g_{st}`$ defined in (33) we have the following implication:

$$\begin{array}{c}((p_1𝐏_{𝐮,𝐉}):[(E_1,E_2_{>0}^1):_{i=1}^nC_i^J=1]):(E_1,E_2_0^1):\hfill \\ \underset{h\mathrm{}}{lim}J=\underset{h\mathrm{}}{lim}g_{st}[E_1,E_2]=\underset{h\mathrm{}}{lim}h\mathrm{cos}c_1f_{st}[\mathrm{tan}c_1]=\underset{h\mathrm{}}{lim}hC_\theta =\mathrm{}.\hfill \end{array}$$
(40)
On the other hand consider a diffusion-limited pathway $`p_2𝐏_{\mathrm{𝐮𝐃},𝐉}`$ where $`J=g[E_1,E_2]`$ as defined by (61). When the system reaches a steady state as defined by (26), we know by (25) that $`\dot{\xi }(t)=0`$ and $`0=I_{in}(t)I_{rev}(t)\dot{o}(t)`$. Since $`\underset{t\mathrm{}}{lim}\dot{o}(t)=J^{}`$, using (28) we have the steady state condition:

$$\begin{array}{c}J^{}=Du(t)Ds_1(t)\text{ and}\hfill \\ J\pm \epsilon =Du(t)Ds_1(t).\hfill \end{array}$$

Since $`\dot{u}(t)=0`$;

$$\begin{array}{c}(p_2𝐏_{\mathrm{𝐮𝐃},𝐉}):(D):(u(0)=c_2):(\epsilon ):\hfill \\ [c_3=Dc_2][J\pm \epsilon =c_3Ds_1(t)].\hfill \end{array}$$
(41)
Recall that $`D0`$ by (28) and $`s_1(t)0`$ by (24). For any finite $`0<\epsilon <\mathrm{}`$ we can define a constant $`c_4`$ such that $`c_4=c_3+\epsilon `$. Consequently

$$Jc_4$$
(42)
where we know that $`0<c_4<\mathrm{}`$. Hence

$$\begin{array}{c}(p_2𝐏_{\mathrm{𝐮𝐃},𝐉}):((E_1,E_2)_{>0}^1:\hfill \\ \underset{h\mathrm{}}{lim}J=\underset{h\mathrm{}}{lim}g[E_1,E_2](c_4<\mathrm{}).\hfill \end{array}$$
(43)
If we now posit a pathway $`p_3`$ that both obeys the flux summation theorem and is diffusion-limited, propositions (40) and (43) lead to the proposition:

$$\begin{array}{c}(p_3𝐏_{\mathrm{𝐮𝐃},𝐉}):\left(\begin{array}{c}[p_3𝐏_{𝐮,𝐉}][(E_1,E_2_0^1):_{i=1}^nC_i^J=1]\hfill \\ [\mathrm{}(c_4<\mathrm{})]\hfill \end{array}\right)\hfill \end{array}$$
(44)
which is a contradiction. Therefore

$$(p𝐏_{\mathrm{𝐮𝐃},𝐉}):[\neg [(E_1,E_2_0^1):\underset{i=1}{\overset{n}{}}C_i^J=1]].$$

$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{14}`$

### B.4 Proof of proposition 16

For a pathway $`p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$, consider a case in which enzyme 2 is saturable such that

$$(E_1_{>0}^1):\left[JE_2k_{cat(2)}J^{}E_2k_{cat(2)}\right],$$
(45)
where for any $`\delta E_1`$, $`J^{}=J+\delta J`$. It follows that if $`E_2`$ is held constant, for any $`\delta E_1`$

$$\frac{\delta J}{\delta E_1}\frac{E_2k_{cat(2)}J}{\delta E_1}.$$
(46)
Multiplying both sides by $`E_1/J`$ we have

$$\frac{E_1}{J}\frac{\delta J}{\delta E_1}\frac{E_1E_2k_{cat(2)}JE_1}{J\delta E_1}.$$
(47)
Hence by (1) we obtain

$$C_1^J\frac{E_1E_2k_{cat(2)}}{J\delta E_1}\frac{E_1}{\delta E_1}.$$
(48)
In a similar fashion if enzyme 1 is saturable and $`E_1`$ held constant, for any $`\delta E_2`$ we have

$$C_2^J\frac{E_1E_2k_{cat(1)}}{J\delta E_2}\frac{E_2}{\delta E_2}.$$
(49)
For each enzyme $`i`$, consider a measure of saturation $`Sat_i`$ such that

$$Sat_i=J/E_ik_{cat(i)}.$$
(50)
It follows that for any enzyme $`i`$,

$$\underset{Sat_i1}{lim}\left(\frac{E_ik_{cat(i)}}{J}\right)=1.$$
(51)
Denoting $`R.H.S.(\text{48})`$ as the right hand side of (48), as enzyme 2 approaches saturation, (51) and (48) imply that

$$\underset{Sat_21}{lim}(R.H.S.(\text{48}))=\frac{E_1}{\delta E_1}\frac{E_1}{\delta E_1}=0.$$
(52)
Similarly, if enzyme 1 is approaching saturation, (51) and (49) imply that

$$\underset{Sat_21}{lim}(R.H.S.(\text{49}))=\frac{E_2}{\delta E_2}\frac{E_2}{\delta E_2}=0.$$
(53)
Consequently when both enzyme 1 and enzyme 2 approach saturation, (52) and (53) imply that

$$(p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}):[(\delta E_1,\delta E_2_{>0}^1):\underset{Sat_1,Sat_21}{lim}\left(\underset{i=1}{\overset{2}{}}C_i^J\right)=0].$$
(54)
$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{16}`$

### B.4 Proof of proposition 17

For a pathway $`p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}`$, consider a case in which enzyme 2 is saturable such that

$$(E_1_{>0}^1):\left[JE_2k_{cat(2)}J^{}E_2k_{cat(2)}\right],$$
(55)
where for any $`\delta E_2`$, $`J^{}=J+\delta J`$.

Given that $`Sat_i=J/E_ik_{cat(i)},`$ if $`E_1`$ is held constant, for any $`\delta E_2<0`$ then

$$\underset{Sat_21}{lim}\left(\frac{E_2k_{cat(2)}}{J}\right)=1\text{ and }\underset{Sat_21}{lim}\left(J\right)=E_2k_{cat(2)}.$$
(56)
Given (55) and (56), for domains in which decreasing $`E_2`$ increases saturation we have:
For all intervals $`a<E_2<b`$ :
If $`\frac{Sat_2}{E_2}0`$ and $`a0`$ and $`\delta E_2<0`$ then:

$$\underset{Sat_21}{lim}\left(\delta J\right)=\delta E_2k_{cat(2)}.$$
(57)
Hence

$$\underset{Sat_21}{lim}\left(C_2^J\right)=\underset{Sat_21}{lim}\left(\frac{\delta J/J}{\delta E/E}\right)=\frac{\delta E_2k_{cat(2)}/k_{cat(2)}}{\delta E_2/E_2}=1.$$
(58)
Similarly if enzyme 1 is saturable we can show that

$$\underset{Sat_11}{lim}\left(C_1^J\right)=1.$$
(59)
Hence for all intervals $`a_i<E_i<b_i`$ :

$$\begin{array}{c}(p𝐏_{𝐮,𝐕_{\mathrm{𝐦𝐚𝐱}},𝐉}):\hfill \\ \left(\begin{array}{c}[i:(\frac{Sat_i}{E_i}0a_i0\delta E_i<0)]\hfill \\ \underset{Sat_1,Sat_21}{lim}\left(C_1^J+C_2^J\right)=2\hfill \end{array}\right)\hfill \end{array}$$
(60)
$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{17}`$

### B.5 Proof of proposition 18

The most straight forward proof of (18) is by direct deduction. To do this we define a set of notations for representing discrete changes in steady state flux due to changes in enzyme concentration.

#### Definition B.4a:

For a given input $`u(t)`$ that is constant, let $`J`$ be the abbreviated notation for the result of a mapping in a function $`g:^2^1`$, with the variables $`E_1`$ and $`E_2`$ as arguments:

$$\begin{array}{c}(p𝐏_{𝐮,𝐉}):\hfill \\ (u(0)^1):((\epsilon >0)_0^1):(g:_0^2^1):\hfill \\ [J=g[E_1,E_2]J=g_\epsilon [u(0),E_1,E_2]]\hfill \end{array}$$
(61)
The discrete operator $`()\stackrel{E_1E_2}{}`$ applied to $`J`$ is defined as:

$$J\stackrel{E_1E_2}{}g[E_1+\delta E_1,E_2+\delta E_2]$$

where the symbol “$``$” denotes identity by definition and for a given function $`g`$, the increments in $`E_1`$ and $`E_2`$ are held constant such that $`\delta E_1=c_1`$, $`\delta E_2=c_2`$ and $`c_1,c_2_0^1`$. Similarly,

$$J\stackrel{E_1}{}g[E_1+\delta E_1,E_2]\text{ and }J\stackrel{E_2}{}g[E_1,E_2+\delta E_2].$$

Note that the operator $`()\stackrel{()}{}`$ can be applied iteratively and is also commutative in its iteration order such that

$$(J\stackrel{E_1}{})\stackrel{E_2}{}=J\stackrel{E_1E_2}{}=(J\stackrel{E_2}{})\stackrel{E_1}{}$$

#### Definition B.5b:

Let the operator $`()\underset{\Delta }{\overset{()}{()}}`$ applied to $`J\stackrel{E_1E_2}{}`$ be defined as

$$J\underset{\Delta }{\overset{E_1E_2}{}}J\stackrel{E_1E_2}{}J$$

where from Definition A1 we get $`(J\stackrel{E_1E_2}{}J)=(g[E_1+\delta E_1,E_2+\delta E_2]g[E_1,E_2])`$. Similarly,

$$J\underset{\Delta }{\overset{E_1}{}}J\stackrel{E_1}{}J\text{ and }J\underset{\Delta }{\overset{E_2}{}}J\stackrel{E_2}{}J.$$

#### Proof Argument:

To prove the relation

$$\begin{array}{c}(E_1,E_2_0^1):\hfill \\ [(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}):\frac{\delta J}{J}=_{i=1}^n\left(\frac{\delta J}{J}\right)_i]\frac{^2J}{E_1E_2}=0,\hfill \end{array}$$
(62)
we start with the left hand side of the relation, which is equation (4) from the Kacser & Burns derivation. Using definitions B.5a and B.5b for a two enzyme pathway, equation (4) is equivalent to

$$\frac{J\underset{\Delta }{\overset{E_1E_2}{}}}{J}=\frac{J\underset{\Delta }{\overset{E_1}{}}}{J}+\frac{J\underset{\Delta }{\overset{E_2}{}}}{J}.$$
(63)
By cancelling the denominator in (63) we proceed with a series of deductions in which the horizontal symbol “$``$” denotes the biconditional relation “ if and only if ” between successive lines. The explanation for each biconditional is given to its right:

$`(E_1,E_2_0^1):\left(\alpha ={\displaystyle \frac{\delta E_1}{E_1}}={\displaystyle \frac{\delta E_2}{E_2}}\right):`$
$`J\underset{\Delta }{\overset{E_1E_2}{}}`$ $`=J\underset{\Delta }{\overset{E_1}{}}+J\underset{\Delta }{\overset{E_2}{}}`$
$`\text{ using  Definitions  B.5a  \&  B.5b}`$
$`g[E_1+\delta E_1,E_2+\delta E_2]g[E_1,E_2]`$ $`=(g[E_1+\delta E_1,E_2]g[E_1,E_2])`$
$`+(g[E_1,E_2+\delta E_2]g[E_1,E_2])`$
$`\text{rearrange}`$
$`(g[E_1+\delta E_1,E_2+\delta E_2]g[E_1+\delta E_1,E_2])`$ $`(g[E_1,E_2+\delta E_2]g[E_1,E_2])=0`$
$`\text{ using  Definitions  B.5a  \&  B.5b}`$
$`(g[E_1,E_2+\delta E_2]g[E_1,E_2])\underset{\Delta }{\overset{E_1}{}}`$ $`=0`$
$`\text{ using  Definitions  B.5a  \&  B.5b}`$
$`(E_1,E_2_0^1):\left(\alpha ={\displaystyle \frac{\delta E_1}{E_1}}={\displaystyle \frac{\delta E_2}{E_2}}\right):`$
$`((J)\underset{\Delta }{\overset{E_2}{}})\underset{\Delta }{\overset{E_1}{}}`$ $`=0.`$
Hence

$$\begin{array}{c}(E_1,E_2_0^1):\left(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}\right):\hfill \\ \left[\frac{\delta J}{J}=_{i=1}^n\left(\frac{\delta J}{J}\right)_i((J)\underset{\Delta }{\overset{E_2}{}})\underset{\Delta }{\overset{E_1}{}}=0\right].\hfill \end{array}$$
(64)
The proposition $`((J)\underset{\Delta }{\overset{E_2}{}})\underset{\Delta }{\overset{E_1}{}}=0`$ is operationally equivalent to the sentence in proposition 14 that proposes the independence of control effects. Next it is simple to show equivalence to a continuous version such that

$$\begin{array}{c}(E_1,E_2_0^1):\hfill \\ [[(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}):((J)\underset{\Delta }{\overset{E_2}{}})\underset{\Delta }{\overset{E_1}{}}=0]\frac{^2J}{E_1E_2}=0].\hfill \end{array}$$
(65)
From the third line of our sequence of biconditional relations we know that the left hand side of (65) has to comply with the equivalence

$$\begin{array}{c}(E_1,E_2_0^1):\left(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}\right):\hfill \\ \left[((J)\underset{\Delta }{\overset{E_2}{}})\underset{\Delta }{\overset{E_1}{}}=0J\stackrel{E_1E_2}{}J\stackrel{E_1}{}J\stackrel{E_2}{}+J=0\right].\hfill \end{array}$$
(66)
We also know by the general definition of partial differentiation that

$$\frac{J}{E_1}=\underset{\delta E_10}{lim}\left(\frac{J\stackrel{E_1}{}J}{\delta E_1}\right)\text{ and }$$

$$\frac{^2J}{E_1E_2}=\underset{\delta E_1,\delta E_20}{lim}\left(\frac{J\stackrel{E_1E_2}{}J\stackrel{E_1}{}J\stackrel{E_2}{}+J}{\delta E_1\delta E_2}\right).$$
(67)
Substituting (66) into (67) we get (65). Using (65) and (64) we get (62).
$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{18}`$

### B.6 Proof of proposition 19

From (11) and (12) we know that the summation theorem implies that

$$\left(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}\right):\frac{\delta J_{1,2}}{J}=\alpha .$$
(68)
Let

$$\delta E_1=\alpha E_1\text{ and }\delta E_2=\alpha E_2.$$
(69)
By definition we know that

$$\begin{array}{c}(E_1,E_2_0^1):(\delta E_1,\delta E_2_{>0}^1):\hfill \\ _{i=1}^2C_i^J=1J=\frac{\delta J_1}{\delta E_1}E_1+\frac{\delta J_2}{\delta E_2}E_2.\hfill \end{array}$$
(70)
Using (68) and (70) we have

$$\delta J_{1,2}=\alpha J=\frac{\delta J_1}{\delta E_1}\alpha E_1+\frac{\delta J_2}{\delta E_2}\alpha E_2.$$
(71)
Using (71) and (69) we have

$$\delta J_{1,2}=\delta J_1+\delta J_2.$$
(72)
By (65) in the proof for proposition 18 we know that

$$\begin{array}{c}(E_1,E_2_0^1):\hfill \\ [[(\alpha =\frac{\delta E_1}{E_1}=\frac{\delta E_2}{E_2}):\delta J_{1,2}=\delta J_1+\delta J_2]\frac{^2J}{E_1E_2}=0].\hfill \end{array}$$
(73)
$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{19}`$

### B.6 Proof of proposition 21

Let $`J=J[E_1,E_2]`$. Given proposition 11 we know that

$$J=hC_\theta \text{ and }J[0,0]=0.$$
(74)
By (74):

$$E_1_0^1:J[E_1,0]=E_1C_{E2=0}.$$
(75)
Let

$$\delta J_{0E1}[E_1,E_2]J[E_1,E_2]J[0,E_2].$$

Using (74) and (75) we have

$$\begin{array}{c}E_1_0^1:\hfill \\ \delta J_{0E1}[E_1,0]=J[E_1,0]J[0,0]=E_1C_{E2=0}.\hfill \end{array}$$
(76)
Separately proposition 19 implies that

$$\begin{array}{c}E_1_0^1:c_\mathrm{\Delta }:E_2_0^1:\hfill \\ \delta J_{0E1}[E_1,E_2]=c_\mathrm{\Delta }.\hfill \end{array}$$
(77)
Using (76) and (77) we have

$$\begin{array}{c}(E_1,E_2_0^1):\hfill \\ \delta J_{0E1}[E_1,E_2]=E_1C_{E2=0}.\hfill \end{array}$$
(78)
Separately using a procedure similar to deriving (75) we have

$$E_2_0^1:J[0,E_2]=E_2C_{E1=0}.$$
(79)
Given

$$J[E_1,E_2]=J[0,E_2]+\delta J_{0E1}[E_1,E_2]$$

and equations (78) and (79) we are left with

$$J[E_1,E_2]=E_2C_{E1=0}+E_1C_{E2=0}.$$
(80)
$`\mathrm{}\mathrm{𝐩𝐫𝐨𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧}\text{21}`$

## C. Derivations for Numerical Analysis

### Derivations for flux surfaces

At steady state the algebraic solution for the state variables of the equations in (22) can be obtained using algebraic solving routines. The solution for the seven internal state variables are in the form of five equations in which there are no explicit solutions for $`e_1`$ and $`e_2`$ such that:

$$\begin{array}{c}s_1=(e_2Q(k_2+k_3)k_5k_7+e_1k_2k_4(Q(k_6+k_7)+e_2k_6k_8))\mathrm{\Theta }\hfill \\ \\ \underset{¯}{es}_1=e_1k_1(e_2Qk_5k_7+e_1k_4(Q(k_6+k_7)+e_2k_6k_8))\mathrm{\Theta }\hfill \\ \\ s_2=e_1k_1k_3(Q(k_6+k_7)+e_2k_6k_8)\mathrm{\Theta }\hfill \\ \\ \underset{¯}{es}_2=e_1e_2k_1k_3k_5(Q+e_2k_8)\mathrm{\Theta }\hfill \\ \\ s_3=e_1e_2k_1k_3k_5k_7\mathrm{\Theta }\hfill \end{array}$$
(81)
where

$$\mathrm{\Theta }=\frac{uD}{e_2Q(k_2+k_3)k_5k_7D+e_1(Qk_2k_4(k_6+k_7)D+e_2(Qk_1k_3k_5k_7+k_2k_4k_6k_8D))}.$$
(82)
In addition recall the stoichiometric constraint that

$$\begin{array}{c}E_1=e_1+\underset{¯}{es}_1\\ E_2=e_2+\underset{¯}{es}_2.\end{array}$$
(83)
Taking the equations (81), (82) and (83) we have a system of seven equations. These can be solved by a convergent Newton-Raphson based numerical algorithm for any case where the set of real-valued parameters

$$\{D,E_1,E_2,Q,k_1,k_2,k_3,k_4,k_5,k_6,k_7,k_8\}$$

is given. The solution for $`\{s_1,s_2,s_3,e_1,\underset{¯}{es}_1,e_2,\underset{¯}{es}_2\}`$ that is arrived at by numerical approximation can then be used to solve for the steady state flux $`J`$. The surfaces in Figure 2 and Figure 7 are numerical solutions of (81), (82) and (83).

The wild-type kinetic parameters used in the numerical procedures were $`k_1=k_5=4\times 10^7M^1sec^1`$ , $`k_2=k_6=4\times 10^2sec^1`$, $`k_3=k_7=7\times 10^2sec^1`$ and $`k_4=k_8=1\times 10^6M^1sec^1.`$ The diffusion constants used are $`D=Q=3\times 10^1sec^1`$. For each enzyme $`K_{eq}=70`$ and for the wild type $`k_{cat}=42,000min^1.`$ Figures 8 and 9 show the steady state solution of all state variables in the pathway $`p_{\text{22}}`$.

For simulations of mutations affecting $`k_{cat}`$ consider enzyme 1 as an example. Kinetic mutations have to occur under the thermodynamic constraint that the equilibrium constant $`K_{eq}=k_1k_3/k_2k_4`$ for the reaction cannot change. Note that $`V_{max}=k_{cat}[E_i]=k_3[E_1]`$. For instance a reduction of $`k_{cat}`$ to 25% of wild type can be modelled by multiplying each of the wild-type rate constants $`k_3`$ and $`k_4`$ by 1/4.

### Failure of a modified summation conjecture

A possible counter to the failure of the summation theorem is to try to include all steps in a pathway, including the non-enzymatic diffusion steps. Hence a conjecture that

$$\underset{i=1}{\overset{n}{}}C_i^J=C_D^J+C_{E1}^J+C_{E1}^J+C_Q^J$$
(84)
Here we show the numerical evidence that such a conjecture does not hold in $`p_{\text{22}}`$. In any case if such a conjecture had been true, it would not have been of much practical value since it would mean that for any pathway deep in metabolism one would have had to take into account every step starting from the initial diffusion barrier. A modification of the proof for proposition (16) can be used to show the failure of the modified summation conjecture in regions of saturation. In Figure 10 we show the numerical failure of the modified summation conjecture for a case where $`\delta E_i=0.02\mu M`$ and $`\delta D/D=\delta Q/Q=0.01`$.


text:
# A Theoretical Search for the Optimum Giant Magnetoresistance

## Abstract

The maximum current-perpendicular-to-plane giant magnetoresistance is searched for in magnetic multilayers made of Co, Ni, and Cu with disorder levels similar to those found in room temperature experiments. The calculation is made possible by a highly optimized linear response code, which uses the impurity averaged Green’s function technique and a 9-band per spin tight binding model. Using simulated annealing, hundreds of different configurations of the atomic layers are examined to find a maximum GMR of 450% in ultrathin Ni/Cu superlattices.

giant magnetoresistance, magnetic multilayers, optimization

The phenomena of the giant magnetoresistance (GMR) has been extensively investigated both theoretically and experimentally since its discovery.Baibich et al. (1988); Binasch et al. (1989) This interest stems from both fundamental questions about magnetism and transport on the nanoscale and also from technological applications like magnetic sensors.Prinz and Hathaway (1995) From this body of work a basic understanding of the physics behind the GMR has emerged. When a magnetic field causes the magnetic domains in a magnetic multilayer or other GMR structure to change their relative orientation, the resistance changes because the bulk and interface scattering rates depend on the electron spin orientation. While this is simple to explain, the actual details are quite complex. Experimentally, it is known that the GMR is very sensitive to atomic composition of the magnetic multilayers. In cases adding a single monolayer can change the magnetoresistance.Parkin (1993) Film growth parameters can similarly change the GMR by an order of magnitude.Colino et al. (1996) The sensitivity to atomic composition and film morphology as well as GMR samples consisting of many atomic layers has made the theory of the GMR challenging (for reviews see Refs. Levy, 1994 and Gijs and Bauer, 1997).

Nonetheless, much progress has been made in elucidating the physics of the GMR in general and in explaining specific experiments. Thus, in this paper we take a different approach. Based on our work and those of other groups, we assume that we have an accurate model of the GMR. Chen et al. (1999) We then ask the question: what is the largest GMR one can obtain for a class of magnetic multilayers with experimentally reasonable disorder parameters? In other words, can one theoretically search the parameter space of a specific class of GMR structures and find the optimum or near optimum atomic configuration? This is a much more challenging problem than computing the GMR for a specific GMR structure because one must perform many individual GMR calculations in the search for the optimum magnetoresistance. Indeed there have been very few calculations like this in condensed matter physics. Franceschetti and Zunger demonstrated the use of simulated annealing to optimize the band gap in semiconductor alloys and superlattices. By searching as much as $`10^4`$ samples, they found the maximum band gap configuration to be very complex.Franceschetti and Zunger (1999) Using a similar approach, Íñiguez and Bellaiche optimized the ground state structural properties and electromachanical response of perovskite alloys.Íñiguez and Bellaiche (2001) This article presents one of the first attempts to optimize a transport property.

Optimizing the transport properties of nanostructures requires a delicate balance between accuracy and speed. To ensure that the optimization results are valid, the model has to include the band structure of the materials, the geometry of the nanostructure, and scattering effects caused by impurities, disorder, and temperature. On the other hand, the calculation of the transport property has to be fast enough that the optimization can be finished in a reasonable amount time. To solve these difficulties, we have developed a highly optimized parallel tight-binding code to compute electron transport in nanostructures using the impurity averaged Green function technique and linear response theory. While the tight binding model we use is not as accurate as other techniques such as density functional theory,Blaas et al. (1999) it does allow us to include realistic band structure efficiently. Similarly, the impurity averaged Green’s function technique allows us to include different kinds of interface and bulk scattering without worrying about the placement of individual defects.

The tight binding Hamiltonian we consider is

$$H=\underset{𝐫,j;𝐫^{},j^{}}{}H_{𝐫,j;𝐫^{},j^{}}c_{𝐫,j}^{}c_{𝐫^{},j},$$
(1)
where $`𝐫`$ and $`𝐫^{}`$ are sites on the lattice. The indices $`j`$ and $`j^{}`$ denote one of the one s, three $`p`$, or five d levels and also the spin orientation. Within a tight binding model, the disorder comes from variations in the matrix elements. In the impurity averaged green function technique, these variations are treated statistically by averaging over an ensemble of different microscopic configurations of the disorder with the same macroscopic properties, e.g., density of impurities. Denoting this average by angular brackets, our unperturbed Hamiltonian is $`H_o=H`$. The deviation from this average for a particular configuration of the disorder is given by $`\delta H=HH_o`$.

Scattering is included through the choice of the self-energy. We use the self-consistent Born approximation with uncorrelated disorder fluctuations. For energy $`\omega `$ the retarded and advanced self-energies, $`\mathrm{\Sigma }^{\text{R}/A}`$, are expressed in terms of the retarded and advanced green functions, $`G^{\text{R}/A}`$ as

$$\mathrm{\Sigma }_{𝐫,j;𝐫,j}^{\text{R}/A}(\omega )=(\delta H_{𝐫,j;𝐫,j})^2G_{𝐫,j;𝐫,j}^{\text{R}/A}(\omega ),$$
(2)
where the green function is in turn expressed in terms of the self-energy via

$$[\omega H_o\mathrm{\Sigma }^{\text{R}/A}(\omega )]G^{\text{R}/A}(\omega )=\mathrm{𝟏}.$$
(3)
This last equation is a matrix equation, where the $`𝐫`$ and $`j`$ indices have been suppressed.

Within the linear response theory, the conductivity is computed via the Kubo formula. For the geometry with the current flowing perpendicular to the planes (CPP) it is necessary to include the vertex correction, $`\mathrm{\Gamma }`$, which is given by the self-consistent equation,

$$𝚪_{𝐫,j;𝐫,j}=(\delta H_{𝐫,j;𝐫,j})^2\left(G^\text{R}(𝐉+𝚪)G^\text{A}\right)|_{𝐫,j;𝐫,j},$$
(4)
where $`𝐉`$ is the current operator,

$$𝐉_{𝐫,j;𝐫^{},j^{}}=i(𝐫𝐫^{})H_{𝐫,j;𝐫^{},j^{}}.$$
(5)
A conserving approximation for the conductivity at $`𝐫`$ is

$`\sigma _{\alpha \beta }(𝐫)`$ $`=`$ $`{\displaystyle \frac{e^2}{h}}\text{Tr}\{J_\alpha (𝐫)G^\text{R}(J_\beta +\mathrm{\Gamma }_\beta )G^\text{A}`$
$``$ $`{\displaystyle \frac{1}{2}}J_\alpha (𝐫)G^\text{R}J_\beta G^\text{R}{\displaystyle \frac{1}{2}}J_\alpha (𝐫)G^\text{A}J_\beta G^\text{A}\},`$
where $`𝐉(𝐫)`$ is the local current density operator at $`𝐫`$. The sum of $`𝐉(𝐫)`$ over position times the lattice unit cell volume is equal to the current operator in Eq. (5). In the CPP geometry this conductivity will be constant because the current density is independent of position.

Letting $`\sigma _\text{F}`$ be the conductivity when the magnetic layers are aligned ferromagnetically and $`\sigma _{\text{A}F}`$ be the conductivity when adjacent magnetic layers are aligned antiferromagnetically, the giant magnetoresistance is defined as $`\text{G}MR=(\sigma _\text{F}\sigma _{\text{A}F})/\sigma _{\text{A}F}`$. Note that in practice the magnetic layers will not necessarily be aligned antiferromagnetically at zero field; however, our simulation does not compute the magnetic domain structure, but only the conductivity. In particular some of the ultrathin magnetic layers which we find with large GMR could very well be ferromagnetically coupled.Parkin et al. (1990) It is possible to achieve antiparallel alignment even with ferromagnetic coupling by going to a spin valve geometry with only two magnetic layers where one of the layers is exchange biased.Dieny et al. (1991)

An optimization study involves many GMR calculations, therefore, the speed for each GMR calculation is crucial. The trace in Eq. (6) is over all positions and atomic levels; however, since the system is periodic, it can be broken up into a sum in k-space and a sum within a unit cell. Thus, the matrices are of order $`9N`$ for each spin, where $`N`$ is the number of atoms in a unit cell. At first glance, the matrix products and inverses in Eqs. (3), (4) and (6) are all $`O(N^3)`$ operations, where $`N`$ is the number of atoms in the unit cell. In addition, both Eqs. (3) and (4) are self-consistent equations. Performing optimization with a self-consistent $`O(N^3)`$ implementation is estimated to take about 5 CPU-years, which is impractical. To overcome this, we have developed an $`O(N)`$ algorithm to solve Eqs. (3) and (6), and an $`O(N^2)`$ algorithm to solve Eq. (4) directly without the need for iteration. For $`N=50`$, this algorithm is more than 100 times faster than the $`O(N^3)`$ algorithm. Running on a Linux cluster of five dual 866MHz Pentium III CPU workstations, it takes about 20 minutes for each GMR calculation and about two weeks for an optimization.

The optimization technique we use is simulated annealing. Since not every configuration qualifies as a GMR structure, the optimization is performed in the constrained subspace consisting of a unit cell with two magnetic layers separated by two non-magnetic spacer layers. After an arbitrarily chosen initial unit cell, each subsequent configuration is generated by one of the following Monte Carlo moves: (i) inserting a monolayer, (ii) removing a monolayer, or (iii) changing the composition of a monolayer. If the GMR of the new configuration is higher than the previous configuration, then it is accepted. On the other hand, if the new GMR is lower, it is only accepted with a probability of $`\mathrm{exp}(\mathrm{\Delta }\text{GMR}/T)`$, where $`\mathrm{\Delta }\text{GMR}`$ is the change in the GMR and $`T`$ is the simulated annealing temperature. The process is continued as the annealing temperature is gradually lowered until there is no further change in the configuration and a final near-global maximum in the GMR is attained.

In the following, we show results of a study that looks for the optimal configuration of superlattices made with Co, Ni, and Cu in conditions similar to GMR experiments at room temperature. To make the study more manageable, we only consider fcc (111) lattices. The minimum thickness of each layer is limited to two monolayers in order to avoid complications such as pin holes. The tight binding parameters in Eq. (1) are obtained from fits to density-functional calculations where up to second-nearest-neighbor hoping energies are included.

Two types of scattering are included: spin-independent scattering and spin-dependent scattering. The spin-independent scattering represents the structural disorder and is included in the calculation by uniformly shifting all the orbital levels at a site: $`(\delta H_{𝐫,j;𝐫,j})^2=0.2\text{e}V^2`$, which corresponds to a Cu resistivity of 3.1 $`\mu \mathrm{\Omega }\text{c}m`$, about twice the value of clean Cu at room temperature. The spin-dependent scattering is assumed to be due to inter-diffusion of atoms near the interface. It is modeled by the shifting of the d-levels only, since the d-levels are the primary difference in the tight binding parameters for these elements. The shift in the d-levels are determined by fits to 4-atom supercell density-functional calculations of an impurity in three host atoms. We then assume an interface mixing percentage of 10% in each of the two monolayers on either side of the interface. The value of $`(\delta H_{𝐫,j;𝐫,j})^2`$ for the majority and minority d-states is equal to 0.29 and 0.84 $`\text{e}V^2`$ at the Co/Cu interface, and 0.03 and 0.11 $`\text{e}V^2`$ at the Ni/Cu interface. Using these parameters, the GMR for a \[Co3/Cu3/Co3/Cu3\] superlattice is 36%, which is in the range of observed CPP-GMR at room temperature.Gijs et al. (1994) Our spin-dependent scattering does not include spin flip scattering, which can be important for short spin flip scattering lengths.Valet and Fert (1993)

The giant magnetoresistance for a simulated annealing run is shown in Fig. 1. The initial configuration has been chosen to be the superlattice \[Co3/Cu3/Co3/Cu3\], with a unit cell made of four layers, each consisting of three monolayers of either Co or Cu. The particular choice of initial configuration does not effect the final result for this calculation. As the annealing temperature $`T`$ decreases from its initial value of $`60\%`$, the GMR becomes larger. The GMR does not increase monotonically because a superlattice that reduces the GMR is accepted with a probability set by the Boltzmann factor, which is crucial for escaping from local maxima in the GMR. As the annealing temperature is lowered, the chance of leaving a local maximum becomes smaller and smaller. At step 89 in this run, the GMR jumps to $`450\%`$ and no other configurations are accepted. The atomic configuration for the highest GMR (d) is \[Ni2/Cu2/Ni2/Cu2\]. This is a very surprising result because most room temperature studies to date have GMR values less than $`100\%`$.Gijs et al. (1994); Gijs and Bauer (1997); Bass and Pratt (1999) From Fig. 1 one can see that the GMR is very sensitive to the atomic configuration. For example, superlattices (b) and (c) have GMR values of 221% and 38%, respectively, yet they differ by only one monolayer from the optimal GMR of superlattice (d), 450%.

To get a better picture of the GMR in configuration space, we have calculated the GMR for 600 randomly generated configurations. The result is shown in Fig. 2 as a histogram. The average GMR is $`22\%`$, and the standard deviation is $`36\%`$. There are two groups of configurations. About $`97\%`$ of the configurations are located near the main peak. They have a GMR less than $`100\%`$. On the other hand, about $`3\%`$ of the configurations have GMR many standard deviations away from the main peak. They scatter in the range from $`100\%`$ to $`450\%`$. Most of the large GMR configurations are ultrathin Ni/Cu superlattices. Only two of the large GMR configurations contain Co. Since the large GMR configurations are rare, they may have been overlooked.

What causes the GMR to be large in some multilayers but not others? The dominant source of scattering in these transition metals is due to s-d hybridization. The d-levels in the majority spin of Ni, Co, and Cu are below the Fermi level; however, there is still a finite d-density of states at the Fermi level due to both s-d hybridization and disorder. In particular disorder can substantially smear the d-density of states, increasing the value at the Fermi level.Chen et al. (1999) Since the d-levels in Ni and Cu are more well matched than those in Co and Cu, the disorder due to interface mixing is smaller in Ni/Cu multilayer than in a Co/Cu multilayer. Thus, the GMR is greater in a Ni/Cu multilayer for the same level of interface mixing because scattering via the s-d hybridization is less.

To see this more clearly, in Fig. 3 we take a representative Ni/Cu sample and a representative Co/Cu sample. The two samples have the same bulk disorder, but the Ni based multilayer actually has more interface mixing (10%) than the one containing Co (5%) to allow the curves to fit clearly on the same plots. As seen in Fig. 3(a), the d-density of states is smaller in the Ni sample. The d-density of states correlates well with the larger conductivity shown in Fig. 3(b) and the larger GMR shown in Fig. 3(c). The fact that there is a simple explanation for this large GMR indicates that it is a robust effect and not strongly dependent on the many assumptions which we have made. Also, it says that for comparable coupling between layers and disorder, the Ni/Cu multilayers will have a larger GMR than the Co/Cu ones.

In this paper we have performed a simulated annealing search for the optimum or maximum giant magnetoresistance by varying the atomic composition of multilayers composed of Ni, Co, and Cu. The disorder parameters in our calculation were chosen to simulate room temperature experimental conditions. After hundreds of GMR calculations for different atomic configurations, a surprisingly large GMR was found for ultrathin Ni/Cu multilayers. The origin of this large GMR of 450% was found to be the fact that the d-levels in Ni and Cu are relatively close in energy. Thus, although we have made a number of approximations in order to perform our calculation in a reasonable amount of time, the large GMR appears to be a robust effect. We believe that this work is part of a large class of problems in tailoring transport properties of nanostructures to optimize an observable such as the magnetoresistance.

###### Acknowledgements.

We would like to thank Fred Sharifi for useful discussions. This work is supported by DOD/AFOSR F49620-1-0026, NSF DMR9357474, and the Research Corp.


text:
# D-Brane Probes of 𝐺₂ Holonomy Manifolds

### Introduction

The concept of the D-brane probe provides the link between gauge theory and geometry. It has proven to be one of the most important ideas among recent developments in string theory, shedding new light on strong coupling gauge theory dualities, providing a physical interpretation of previously obscure algebraic-geometric constructions, and resulting in examples of the AdS/CFT correspondence with increasingly realistic gauge and matter content.

The purpose of this letter is to describe some results extending this enterprise to D2-brane probes of non-compact manifolds of $`G_2`$ holonomy. These spaces are of interest for both theoretical and phenomenological reasons bob ; amv ; aw ; bobed ; decon . For orbifold backgrounds, and their partial resolutions, there are well known procedures for determining the gauge theory living on the probe world-volume dm ; drm . (See tsimpis  for examples of $`G_2`$ holonomy orbifolds). However, recent attention has focused on asymptotically conical $`G_2`$ manifolds for which the orbifold prescription is not relevant. Here we suggest a different approach to understand aspects of the probe theory, using mirror symmetry of three dimensional gauge theories is .

Recall that mirror symmetry is a duality between a pair of three-dimensional field theories which, among other things, interchanges the Coulomb and Higgs branches. The basic observation is that mirror symmetry on an M2-brane probe has a simple interpretation as an “M-theory flip” pz :

$$\begin{array}{ccccc}& & \text{M-theory on }X& & \\ & & & & \\ \text{IIA with}& & & & \\ \text{D6-branes}& & & & \text{IIA on }X\\ & & & & \\ \text{Coulomb}& & \text{mirror symmetry}& & \text{Higgs}\\ \text{branch}& & & & \text{branch}\end{array}$$

In order to explain this duality, let us consider a membrane probe on a (singular) space $`X`$ with $`G_2`$ holonomy. We can perform a reduction from M-theory to IIA string theory in two different ways. First, we can reduce on a circle transverse to both the manifold $`X`$, and to the membrane. This takes us to IIA theory on $`X`$. The space $`X`$ is reproduced as a Higgs branch in the D2-brane world-volume theory, in a way reminiscent to the hyperKähler quotient construction. By analogy, we call it the $`G_2`$ quotient construction.

On the other hand, we can reduce on an $`𝐒^1`$ contained within $`X`$. There are many ways to choose the $`𝐒^1`$, which lead to different IIA backgrounds with D6-branes and/or Ramond-Ramond fluxes. An illustrative example was described in amv . One might ask: “What is the natural choice of the M-theory circle?”. One simple choice is to require $`X/U(1)^6`$. If such a quotient exists it gives, after reduction to type IIA theory, a configuration of D6-branes in a (topologically) flat space-time. The positions of the D6-branes correspond to the fixed points of the circle action. In contrast to the previous reduction, the geometry of $`X`$ is entirely encoded in the configuration of D6-branes, rather than in the geometry of space-time aw ; GS . The M-theory geometry $`X`$ is then re-constructed as the Coulomb branch of the world-volume theory on the D2-brane. When the space $`X`$ develops a conical singularity, the configuration of D6-branes also becomes singular. In particular, in some cases of interest it degenerates into a collection of flat D6-branes intersecting at the suitable angles BDL . For such models, it becomes a simple exercise to derive the three-dimensional gauge theory on the D2-brane probe. Then, using mirror symmetry, one can obtain the theory on the probe of $`X`$.

The simplest illustration of this method was given in pz , where a single D2-brane probe of $`N`$ parallel D6-branes was used to re-derive the original $`𝒩=4`$ mirror pairs of Intriligator and Seiberg is . Note, however, that we here we run the logic of pz  in reverse: we use mirror symmetry, derived through independent techniques, to derive the gauge theory on the probe of $`X`$.

This letter contains only the barest details of our method. Many further results and examples, including the derivation of the mirror pairs and the interesting subtleties involved, as well as applications to $`SU(3)`$, $`Sp(2)`$ and $`Spin(7)`$ holonomy manifolds can be found in gt .

### Mirror Symmetry

Let us firstly describe the mirror pairs that will be our tool in understanding the probe theories. Our mirror pairs have $`𝒩=1`$ supersymmetry, and are derived in gt  from deforming $`𝒩=4`$ mirror pairs using both field theory and string theory techniques. Of course with such little supersymmetry ($`𝒩=1`$ means 2 supercharges) we have little control over the strong coupling dynamics and must be wary of any conjectured duality. Our only savior is parity symmetry which may be used to prohibit the lifting of certain vacuum moduli spaces ahw ; GK . We hope that the success of our mirror pairs in describing manifolds of $`G_2`$ holonomy goes some way towards convincing the reader of their utility.

The mirror pairs preserve only $`𝒩=1`$ supersymmetry, and are given by

$`\mathrm{𝐓𝐡𝐞𝐨𝐫𝐲}𝐀:`$ $`U(1)^r\text{with }k\text{ scalars and }N\text{ hypers}`$
$`\mathrm{𝐓𝐡𝐞𝐨𝐫𝐲}𝐁:`$ $`U(1)^{Nr}\text{with }(3Nk)\text{ scalars and}`$
$`N`$ hypermultiplets
The abelian vector multiplets contain only a photon and a Majorana spinor, while the scalar multiplets, which we shall denote as $`\mathrm{\Phi }`$, contain a single real scalar and a Majorana fermion. In contrast, the hypermultiplets fill out representations of the $`𝒩=4`$ algebra: they each contain four Majorana fermions and two complex scalars, $`q`$ and $`\stackrel{~}{q}`$. We write the superfield as a doublet, $`W=(Q,\stackrel{~}{Q}^{})^T`$. The chiral multiplets $`Q`$ and $`\stackrel{~}{Q}`$ carry conjugate charges under the gauge group. For Theory A, we denote the charge of the hypermultiplets as $`R_i^a`$, while for Theory B it is $`\widehat{R}_i^p`$; $`i=1,\mathrm{},N;a=1,\mathrm{},r;p=1,\mathrm{},Nr`$. Each of these matrices are assumed to be of maximal rank, and mirror symmetry requires

$$\underset{i=1}{\overset{N}{}}R_i^a\widehat{R}_i^p=0a,p$$
(1)
In $`𝒩=1`$ theories, there are no holomorphic luxuries and interactions are determined in terms of a real superpotential, $`f`$. For Theory A, this superpotential has the cubic form associated with the $`𝒩=4`$ theories, and is determined by a triplet of $`k\times N`$ matrices $`T_c`$, $`c=1,2,3`$

$$f=\underset{i=1}{\overset{N}{}}\underset{c=1}{\overset{3}{}}\underset{\alpha =1}{\overset{k}{}}W_i^{}\tau ^cW_iT_{c,i}^\alpha \mathrm{\Phi }_\alpha $$
(2)
where $`\alpha =1,\mathrm{},k`$ and $`\tau ^c`$ are the three Pauli matrices. A similar coupling exists for Theory B, now with the triplet of $`(3Nk)\times N`$ coupling matrices $`\widehat{T}_c`$, satisfying

$$\underset{i=1}{\overset{N}{}}\underset{c=1}{\overset{3}{}}T_{c,i}^\alpha \widehat{T}_{c,i}^\rho =0\alpha ,\rho $$
(3)
Further details of these theories, together with the methods used to derive them, can be found in gt . Here let us restrict ourselves to a few comments. The Coulomb branch of Theory A has dimension $`(N+r)`$, which coincides with the dimension of the Higgs branch of Theory B. (The converse also holds). Mass and FI parameters may be added to both theories, partially lifting some branches of vacua, and the mirror map for these deformations is known.

### The $`G_2`$ Quotient Construction

Let us now apply the mirror pairs described above to a D2-brane probe of a D6-brane background. We take $`i=1,\mathrm{},N`$, flat D6-branes, each of which has spatial world-volume direction,

$`D6_i:`$ $`123[47]_{\theta _1^i}[58]_{\theta _2^i}[69]_{\theta _3^i}`$
The D6-branes lie on a special lagrangian locus if each rotation is contained in $`SU(3)`$ BDL  or, more simply, if

$$\theta _1^i\pm \theta _2^i\pm \theta _3^i=0\mathrm{mod}2\pi i$$
(4)
ensuring that $`𝒩=1`$ supersymmetry (4 supercharges) is preserved on their common world-volume. (For non-generic angles, more supersymmetry may be preserved. We will assume this is not the case).

As described in the introduction, we probe this configuration with a D2-brane lying in the $`x^1x^2`$ plane. This breaks supersymmetry by a further half, resulting in a $`d=2+1`$ dimensional world-volume theory with $`𝒩=1`$ supersymmetry (2 supercharges). For the singular case of intersecting, flat D6-branes, the theory on the D2-brane probe is simple to write down. The 2-2 strings give rise to the usual gauge field and seven scalars. Of these, there is one free $`𝒩=1`$ scalar multiplet parameterizing motion in the $`x^3`$ direction common to all D6-branes. Further fields arise from the 2-6 strings. These give rise to $`N`$ hypermultiplets. Thus, we have the interacting $`𝒩=1`$ supersymmetric theory on the probe,

Theory A: $`U(1)`$ with 6 scalar multiplets and N hypermultiplets

where each hypermultiplet has charge $`+1`$ under the gauge field. The couplings of the hypermultiplets to the scalar multiplets are determined by the geometry of the D6-branes: each hypermultiplet couples minimally to the three scalar fields orthogonal to the corresponding D6-brane. If we define the scalar fields $`\varphi _\alpha =x^{\alpha +3}`$, $`\alpha =1,\mathrm{},6`$, then the superpotential is of the form (2) with the couplings determined by the D6-brane orientations,

$$T_{c,i}^\alpha =\mathrm{sin}\theta _c^i\varphi _\alpha \delta _{c,\alpha }+\mathrm{cos}\theta _c^i\varphi _\alpha \delta _{c,\alpha 3}c=1,2,3$$
(5)
From the IIA space-time picture, we are lead to the natural conjecture that the Coulomb branch of this theory, parameterized by the six real scalars $`\varphi _\alpha `$, together with the dual photon $`\sigma `$, is a seven dimensional manifold $`X`$ that admits a metric of $`G_2`$ holonomy. However, this description of $`X`$ in terms of Coulomb branch variables is not overly useful. In particular, the isometries of $`X`$ are lost in the reduction to IIA, and are only expected to be recovered as isometries of the Coulomb branch in the strong coupling limit. It would be desirable to have an algebraic description of $`X`$, in which the symmetries are manifest. This is exactly what the mirror theory provides for us.

Since Theory A is in the class of theories discussed above, we may simply write down the mirror theory whose Higgs branch is conjectured to give the $`G_2`$ manifold $`X`$,

Theory B: $`U(1)^{N1}`$ with $`3(N2)`$ scalar and $`N`$ hypermultiplets

The gauge couplings are determined by the $`A_{N1}`$ quiver diagram: i.e the $`i^{\mathrm{th}}`$ gauge group acts on the $`i^{\mathrm{th}}`$ hypermultiplet with charge $`+1`$, and the $`(i+1)^{\mathrm{th}}`$ hypermultiplet with charge $`1`$. All other hypermultiplets are neutral. The Yukawa terms are of the form (2), with the triplet of coupling matrices $`\widehat{T}`$ determined by (3). The Higgs branch of this theory is parameterized by $`w_i=(q_i,\stackrel{~}{q}_i^{})^T`$, the $`N`$ doublets of complex scalars in the hypermultiplets. These are constrained by the $`3(N2)`$ D-terms, modulo $`(N1)`$ $`U(1)`$ gauge quotients,

$$\underset{i,c}{}\widehat{T}_{c,i}^\rho w_i^{}\tau ^cw_i=0\rho =1,\mathrm{},3(N2)$$
(6)
This quotient construction yields a conical manifold, which is expected to admit a metric of $`G_2`$ holonomy. In some cases the conical singularity may be (partially) resolved by adding constants to the right-hand side of (6). This blows up two-cycles and, in the IIA picture, corresponds to translating the D6-branes in the $`x^4x^9`$ directions. Note that when the Yukawa matrices $`\widehat{T}`$ fall into suitable $`SU(2)`$ triplets, the above method coincides with the toric hyperKähler quotient construction, supplemented by a further quotient by a tri-holomorphic isometry to yield a manifold of dimension seven. This is the construction discussed by Acharya and Witten bobed . However, in general, the charges in (6) differ.

The above theory may also be considered as a $`𝒩=(1,1)`$ supersymmetric linear sigma model in $`d=1+1`$ dimensions, cf. av . However, in the absence of something akin to Yau’s theorem, we cannot be sure that the Ricci flat metric to which the theory flows has $`G_2`$ holonomy.

### An Example

Let us now examine the $`G_2`$-quotient construction applied to a specific example. Our choice for consideration is the $`G_2`$ manifold $`X`$ given by the cone over the flag manifold $`SU(3)/U(1)^2`$ BS ; gary . This example was also discussed in detail by Atiyah and Witten aw . They show that, with a suitable choice of M-theory circle, $`X`$ can be reduced to three, flat, intersecting D6-branes in type IIA string theory. The symmetry of $`X`$ (to be discussed below), together with the special lagrangian condition (4) determines the angles of these three branes to be $`\theta _c^1=0`$, $`\theta _c^2=2\pi /3`$ and $`\theta _c^3=4\pi /3`$ for each $`c`$. The configuration is drawn in Figure 1.

In order to make the symmetries of the configuration manifest, we define two triplets of scalars, $`\stackrel{}{\varphi }_1=(x^7,x^8,x^9)^T`$ and $`\stackrel{}{\varphi }_2=(x^4,x^5,x^6)^T`$ in terms of which the orientation of the $`i^{\mathrm{th}}`$ D6-brane can be described by the set of linear equations.

$`D6_1:`$ $`\stackrel{}{\varphi }_1=0`$
$`D6_2:`$ $`\frac{1}{2}\stackrel{}{\varphi }_1+\frac{\sqrt{3}}{2}\stackrel{}{\varphi }_2=0`$
$`D6_3:`$ $`\frac{1}{2}\stackrel{}{\varphi }_1+\frac{\sqrt{3}}{2}\stackrel{}{\varphi }_2=0`$
The original $`G_2`$ holonomy manifold $`X`$ enjoyed an $`SU(3)`$ continuous isometry. It’s not surprising that, upon taking the quotient to IIA string theory, this isometry is partially lost. In fact, the D6-brane background has only a $`SU(2)`$ symmetry, under which each $`\stackrel{}{\varphi }_a`$ transforms as a triplet. The M-theory circle itself provides one further, hidden, $`U(1)`$ action. We therefore conclude that the reduction to IIA string theory has broken the isometry group to $`SU(3)SU(2)\times U(1)`$. Now, let us introduce a probe D2-brane in this background, and look at the $`𝒩=1`$ gauge theory on its world-volume.

Theory A: $`U(1)`$ with 6 scalars and 3 hypermultiplets

As described above, the 6 scalar multiplets combine into two triplets whose interactions with the hypermultiplets are of the form (2), where the interaction matrices are determined by (5). The Coulomb branch of this theory is parameterized by the six scalars, together with the dual photon. It has the $`SU(2)\times U(1)`$ isometry group, which is expected to be enhanced to the full $`SU(3)`$ only in the strong coupling, infra-red limit.

Using the results described earlier, the mirror theory is the $`𝒩=1`$ gauge theory with matter content,

Theory B: $`U(1)^2`$ with 3 scalar and 3 hypermultiplets

The charges of the three hypermultiplets under the $`U(1)^2`$ gauge group are $`(+1,1,0)`$ and $`(0,+1,1)`$. The three scalars couple through the usual superpotential (2), with interactions determined by (3) and (5) to be $`\widehat{T}_{c,i}^\rho =\delta _c^\rho `$ for all $`i`$. Let us examine the Higgs branch of this theory. The superpotential provides 3 real constraints on the 12 real scalar fields contained in the hypermultiplets. After dividing by the gauge group, we are left with a Higgs branch of dimension 7, as required. The constraints are,

$`{\displaystyle \underset{i=1}{\overset{3}{}}}|q_i|^2|\stackrel{~}{q}_i|^2=0,{\displaystyle \underset{i=1}{\overset{3}{}}}\stackrel{~}{q}_iq_i=0`$ (7)
Firstly notice that this space has a manifest $`SU(3)`$ symmetry, thus recovering the full isometry group of $`X`$. It is not difficult to further show that the space is indeed isomorphic to the cone over $`SU(3)`$, ensuring that the full Higgs branch is the flag manifold $`SU(3)/U(1)^2`$.

There is a single normalizable deformation of this space, which yields a smooth $`G_2`$ manifold:

$$X^3\times 𝐏^\mathrm{𝟐}$$

In the D6-brane picture, the singularity is resolved by deforming the singular locus of flat, intersecting D6-branes into a smooth special Lagrangian curve $`L^3`$:

$$L\times 𝐒^\mathrm{𝟐}^\mathrm{𝟑}$$

In the present case this deformation involves only two out of the three D6-branes. To see this more explicitly, let us choose the first and second D6-branes, which deform to lie on the special lagrangian curve:

$$\stackrel{}{\varphi }_1\stackrel{}{\varphi }_2=|\stackrel{}{\varphi }_1||\stackrel{}{\varphi }_2|,|\stackrel{}{\varphi }_1|(3|\stackrel{}{\varphi }_1|^2|\stackrel{}{\varphi }_2|^2)=\rho $$
(8)
This curve has a remarkable property: it creates a hole through which the remaining D6-brane can pass, see Figure 1. Therefore, it suffices to deform only two of the three D6-branes in order to completely remove the conical singularity. Of course, one has three different ways to pick a pair of D6-branes, leading to three different resolutions of the space, meeting at a singular point. This is precisely the picture suggested in aw .

It is natural to ask how the probe theory responds to such a deformation. From the perspective of Theory A, one can show that there is essentially a unique deformation consistent with all the symmetries of the model; it is a Yukawa term coupling a pair of hypermultiplets. Moreover, the locus of zeroes of the fermion mass matrix has the same topology as the locus (8). For more details, see gt .

## Acknowledgements

We are grateful to B. Acharya, M. Aganagic, N. Constable, A. Hanany, J. Sparks, N. Seiberg, M. Strassler, C. Vafa and E. Witten for useful discussions. This research was conducted during the period S.G. served as a Clay Mathematics Institute Long-Term Prize Fellow. The work of S.G. is also supported in part by grant RFBR No. 01-02-17488, and the Russian President’s grant No. 00-15-99296. D.T. is a supported by a Pappalardo fellowship and, in part, by funds provided by the U.S. Department of Energy (D.O.E.) under cooperative research agreement #DF-FC02-94ER40818.