Effects of epistasis on phenotypic robustness in metabolic pathways

Mathematical Biosciences 184 (2003) 27–51 www.elsevier.com/locate/mbs

Effects of epistasis on phenotypic robustness in metabolic pathways Homayoun Bagheri-Chaichian

a,b,*

, Joachim Hermisson b, Juozas R. Vaisnys G€ unter P. Wagner b

b,c

,

a

b

Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA Department of Ecology and Evolutionary Biology, Yale University, New Haven, CT 06520-8106, USA c Department of Electrical Engineering, Yale University, New Haven, CT 06520, USA Received 6 September 2002; received in revised form 11 February 2003; accepted 18 March 2003

Abstract It is an open question whether phenomena such as phenotypic robustness to mutation evolve as adaptations or are simply an inherent property of genetic systems. As a case study, we examine this question with regard to dominance in metabolic physiology. Traditionally the conclusion that has been derived from Metabolic Control Analysis has been that dominance is an inevitable property of multi-enzyme systems and hence does not require an evolutionary explanation. This view is based on a mathematical result commonly referred to as the flux summation theorem. However it is shown here that for mutations involving finite changes (of any magnitude) in enzyme concentration, the flux summation theorem can only hold in a very restricted set of conditions. Using both analytical and simulation results we show that for finite changes, the summation theorem is only valid in cases where the relationship between genotype and phenotype is linear and devoid of non-linearities in the form of epistasis. Such an absence of epistasis is unlikely in metabolic systems. As an example, we show that epistasis can arise in scenarios where we assume generic non-linearities such as those caused by enzyme saturation. In such cases dominance levels can be modified by mutations that affect saturation levels. The implication is that dominance is not a necessary property of metabolic systems and that it can be subject to evolutionary modification. Ó 2003 Elsevier Science Inc. All rights reserved. Keywords: Robustness; Dominance; Metabolic control analysis; Epistasis; Summation theorem

*

Corresponding author. Address: Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA. Tel.: +1-505 984 8800; fax: +1-505 982 0565. E-mail address: [email protected] (H. Bagheri-Chaichian). 0025-5564/03/$ - see front matter Ó 2003 Elsevier Science Inc. All rights reserved. doi:10.1016/S0025-5564(03)00057-9

28

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

1. Introduction A recurrent debate in theoretical population genetics has been the nature of dominance and its evolution [1]. The main question has centered around whether the widely observed dominance of wild-type phenotypes to their mutant counterparts is an accidental property of genetic systems or whether it is an evolved adaptation [2]. This is an important question because it addresses whether the effects of mutations could be altered by natural selection in such a way that mutations would not be detrimental to individuals in a population. In the past, conclusions derived on the issue of dominance have been mixed. There are two parts to the question of dominance evolution, and they have been usually addressed in nonoverlapping disciplines with different intellectual histories. One question is whether dominance properties of a genetic system could be altered by mutations. The second is whether population conditions exist for the selection of dominance, since these conditions do not necessarily follow even if dominance properties could be altered by mutations. As far as dominance modification was concerned, many early experiments confirmed that this was possible in the laboratory [3–7]. However population genetics models indicated that the selection coefficient for dominance modifiers would not be strong enough to overcome drift in populations with normal mutation rates [8–11]. Hence the genetic studies indicated that dominance could be modified by mutation while the population studies indicated that it could not be selected for. The only exception were cases where the frequency of heterozygotes were to be high in a population due to a balanced polymorphism. Under such circumstances the consensus was that the selection coefficients were high enough for the evolution of dominance to occur [11–16]. The 1980s witnessed a further influx of dissenting views on this topic. Based on their studies of multi-enzyme systems, theoreticians in the field of Metabolic Control Analysis (MCA) concluded that in the realm of metabolic physiology, dominance was a necessary property of a pathway and could not be significantly altered by mutation [17]. Conversely in population genetics, new models indicated that dominance could be selected in populations that started far from equilibrium and had a high initial frequency of deleterious mutants [18–20]. The result was a divergence into two opposing positions. One position based on the convergence of MCA and results from equilibrium population genetics, held that dominance could not evolve [21–26]. Another, based on several approaches, such as observations on dominance modification [1,27,28], results from non-equilibrium population genetics [18–20], expected properties of non-linear systems [29–32], and examples involving balanced polymorphisms [16,33] was that dominance could evolve. In order to understand dominance we have to understand how phenotypic robustness to mutation can arise. Robustness is a term used for the class of phenomena whereby a property (such as the phenotype) is not sensitive to perturbations in an underlying variable (for example, genetic variation). As such, phenomena such as genetic canalization, dominance and redundancy are specific examples of phenotypic robustness to mutations [34–36]. In this paper we will examine the robustness properties of pathway flux (as the phenotype) with respect to mutations that change the underlying enzyme concentrations. In this context, the question of dominance evolution is in essence a question on whether phenotypic robustness to mutations could result from selection for such robustness [35,37]. As an illustration, consider a genetic locus ÔAÕ, that codes for a metabolic enzyme. Let ÔAÕ designate the wild type allele and ÔaÕ the null mutant allele. Fur-

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

29

thermore consider a case in which the concentration of functional enzyme in an AA homozygote is higher than in an Aa heterozygote; whose concentration is in turn higher than an aa homozygote (which has zero concentration). If metabolic flux for the AA homozygote is robust with respect to decreases in enzyme concentration, then the flux in an Aa heterozygote will not be too different than that of the AA homozygote. Under such circumstances the wild type phenotype associated with the AA genotype is deemed to be dominant with respect to the mutant phenotype associated with the aa genotype. Hence dominance of the wild type could evolve if robustness could evolve. Nonetheless it is imperative to note that the form of dominance illustrated here is simply the mechanistic model that has been debated most within the topic of dominance evolution [1,8–10,17,24–26,29]. Dominance can arise through other mechanisms, an example being compensatory feedback regulation [31]. In recent years there has been a resurgence of interest in the issue of phenotypic robustness to mutations [38–46]. This renewed interest is partly due to the fact that genetic experiments are indicating that many knock-out mutations have little effect on phenotype [42] and that developmental systems show a high degree of stability with respect to perturbations [41,43]. This renewed interest brings back to the foreground the necessity of resolving the problem of dominance evolution. In this paper we will address one aspect of the problem: the claim in MCA that dominance is an inevitable property of multi-enzyme systems. We will use the same theoretical apparatus used in MCA to show that in fact robustness properties in metabolic systems are quite malleable and not an inevitable expectation.

2. Revisiting the MCA argument for the inevitable expectation of dominance Elucidating how changes in enzyme properties alter the physiological phenotype is a key component for understanding the relation between genotype and phenotype. The two main approaches to this problem have developed into fields known as Metabolic Control Analysis (MCA) [21,47–49] and Biochemical Systems Theory (BST) [50,51]. 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 with respect to phenotypic robustness 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. First, the summation theorem implies that in general the control of flux in a pathway is shared between enzymes. Hence rate limiting enzymes are rare [22]. Secondly, the theorem implies that on the average, mutations that change enzyme concentrations will have a small effect on flux. If true, this implication would have important consequences in evolutionary theory and genetics. Many phenotypes are dependent on metabolism. For such phenotypes, phenomena that fall under the rubric of phenotypic robustness to mutations (such as selective neutrality, canalization and dominance) would be an inevitable outcome of the underlying metabolism and not a result of evolution. Kacser and Burns [17,24,25] made such an argument with regard to the case of dominance, in particular stating that dominance is an inevitable property of metabolic pathways and not a result of evolution.

30

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

At the core of the MCA approach is a measure of the control exerted by each enzyme on flux through a pathway. The flux control coefficient CiJ measures how important each enzyme is in its ability to affect steady-state flux via changes in enzyme concentration. In its original formulation [47] the control coefficient was defined as CiJ ¼

dJi =J ; dEi =Ei

ð1Þ

where J is the steady-state flux of metabolites through the pathway (net rate of product formation at the end of a pathway), Ei is the concentration of enzyme i, dEi is a finite change in concentration of enzyme i and dJi is the resultant change in flux. Hence CiJ is a non-dimensionalized sensitivity measure that shows the ratio of proportional change in flux to proportional change in enzyme concentration. In many models enzyme activity is used instead of concentration. The two uses are equivalent as long as enzyme activity is proportional to enzyme concentration. The flux summation theorem states that the sum of the control coefficients in a pathway with n enzymes equals one: n X

CiJ ¼ 1:

ð2Þ

i¼1

Eq. (2) implies that the average expectation for most enzymes will be that their control coefficient will be on the order of 1=n (though of course it can be much higher for some). As n increases, the control coefficient for most enzymes will get smaller on average. Hence if (2) were to be true most enzymes would have small effects on flux and robustness would be an inevitable property of metabolic pathways. Kacser and Burns connected this robustness argument to the topic of dominance. They argued that in mutant heterozygotes the concentration of functional enzyme for a given locus will be lower than in the the wild type homozygotes (since half the transcripts will be dysfunctional). However if a pathway exhibits phenotypic robustness to enzyme concentrations then the phenotype of the heterozygote would not be very different from the wild type homozygote. A different summation relation can be formulated in terms of a more general definition of control coefficients in which control is defined in terms of local reaction rates and underlying parameters such that CiJ ¼

dJi =J ; dvi =vi

ð3Þ

where vi is the local rate in isolation of the reaction catalyzed by enzyme i [52,53]. However in the case of dominance all arguments have been made in terms of enzyme activity and use equations of the form portrayed in (1) and (2) [17,24,25,54–56]. Similarly, experimental observations of phenotypic robustness have been interpreted in terms of Kacser and BurnÕs dominance argument based on enzyme activities [23,57–61]. Nonetheless it is worthwhile noting that for the common Michaelis–Menten case – in which if the substrate concentrations are held fixed then the isolated reaction rates are linearly dependent on enzyme concentrations – then the two formulations are equivalent. Note that if vi ¼ cEi , where c is a constant, then dvi =vi ¼ dEi =Ei and hence (3) reduces to (1) [52,53,62]. In cases where vi 6¼ cEi for fixed substrate concentrations then a different summation relation would have to be used [52,63,64].

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

31

Theoretical results that contradict conclusions derived within MCA with regard to dominance have been posited several times [29,30,51,65]. First, Cornish-Bowden [65] showed that in a sequential pathway if the maximal rate Vmax 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 [66]. A second set of objections was then put forth by Savageau and Sorribas [29,51]. 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 some MCA proponents [67]. The basis of the rejection is unclear and results in more recent work are consistent with this objection [30,31,68]. It is important to note that most scientists use the continuous version of (1). However for the case of phenotypic robustness the finite version is the relevant form. Practitioners in MCA have long recognized problems of dealing with finite changes, or what is sometimes referred to as Ôlarge changesÕ [22,53,62,69–71]. However these have never been re-oriented with respect to the problem of dominance and phenotypic robustness. Mutations necessarily involve finite changes in enzyme activity. Hence, whether a continuous version of the summation theorem holds or not (which it would in cases where the flux function is homogeneous, see [72]) does not say anything about whether a pathway will be robust or exhibit dominance with respect to finite changes of any magnitude. In order for a pathway to be robust with respect to finite changes, the finite version of the summation theorem has to hold. In this paper we will address the shortcomings of the summation theorem with respect to the issue of phenotypic robustness. We prove that the summation theorem would hold for finite changes in enzyme concentration if and only if flux is a linear function of enzyme concentrations. This implies that if a multi-enzyme system allows for non-linear effects such as epistasis then the finite version of the summation theorem does not hold. For example if we allow for the possibility of enzyme saturation then the sum of control coefficients can be as high as the the number of enzymes in a pathway. Hence there are no a priori constraints that would require the magnitude of mutational effects to be small. The implication is that phenotypic robustness is not an inevitable property of multi-enzyme systems and is instead an evolvable property. 3. Properties of two enzyme pathways Consider the class of sequential two enzyme pathways that transform an input substrate s1 to a final product s3 which is irreversibly removed by a linear sink step to an output o: Enzyme 1

Enzyme 2

sink

s1 ! s2 ! s3 ! o: MCA applies to the sub-class of pathways whereupon if the input substrate s1 is kept constant, the pathway can reach a steady-state flux J such that J ¼ do=dt and dsi =dt ¼ 0 for all i 2 f1; 2; 3g. Furthermore it is assumed that for every input s1 we can construct a function g : R2 7! R1 such that ð4Þ J ¼ gðE1 ; E2 Þ; where E1 and E2 are the total concentrations of enzymes 1 and 2 respectively.

32

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

As one of the simplest cases of a multi-enzyme system, we examine some of the properties of two enzyme pathways whose flux can be determined by a function of the form g. We shall mainly concentrate on the relation of g with respect to the sum of control coefficients for finite changes in enzyme concentration. We list these properties as a series of propositions whose proofs are given in Appendix A. 3.1. The effects of epistasis At its core, the question of dominance modification or robustness modification is a question about epistasis [31,32,37,73,74]. We are asking whether alleles at one locus can alter the phenotypic effect of alleles at a different locus. In a similar fashion 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. Epistasis is the phenotypic manifestation of interactions and can be measured as the mutual effect that two enzymes have on each-otherÕs ability to affect pathway flux. The condition of no epistasis can be defined as the situation in which o2 J =oE1 oE2 ¼ 0. One can prove that for a two enzyme pathway, the summation theorem implies that the change in flux due to a change of enzyme concentration at one locus is independent of enzyme concentration at the second locus in the pathway. In other words for any given set of enzyme concentrations, if C1J þ C2J ¼ 1 for all ðdE1 ; dE2 Þ then

o2 J ¼ 0: oE1 oE2

ð5Þ

Hence the summation theorem implies the inexistence of epistasis. It immediately follows that o2 J 6¼ 0; oE1 oE2 there is at least one pair ðdE1 ; dE2 Þ such that C1J þ C2J 6¼ 1: for any ðE1 ; E2 Þ where

ð6Þ

In any regime where epistasis is present the finite version of the summation theorem is violated. Hence whether dominance could evolve or not in a metabolic pathway revolves around whether metabolism exhibits epistasis. If there is no epistasis then dominance effects cannot be modified. On the other hand if epistasis is present, then dominance can be modified.

3.2. Summation to unity for finite changes implies a linear flux function In addition to the inexistence of epistasis, one can further show that for any enzyme i, if C1J þ C2J ¼ 1 for all ðdE1 ; dE2 Þ then

o2 J ¼ 0: oEi2

ð7Þ

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

33

Eq. (7) combined with (5) implies that the finite summation theorem holds if and only if flux is a linear function of enzyme concentrations. In other words C1J þ C2J ¼ 1 for all ðdE1 ; dE2 Þ if and only if J ¼ E1 c1 þ E2 c2 ;

ð8Þ

where c1 and c2 are constants. This is a simple but compelling result. The claims that dominance could not evolve would only hold if metabolism could be represented by a linear function. This is very unlikely.

3.3. The effects of non-linearity as a result of enzyme saturation What are the likely causes of non-linearity and epistasis in metabolic pathways? One obvious candidate is the capability of enzymes to be saturated. The possibility of saturation is a very general property of enzyme catalyzed reactions. In terms of its effects on flux in a sequential pathway, it means that for every saturable enzyme i, there exists a constant kcatðiÞ (pertaining to the rate limiting step of catalysis) such that for any concentration of Ei : J 6 Ei kcatðiÞ :

ð9Þ

Accordingly, for each enzyme i we can consider a measure of saturation Sati such that Sati ¼ J =Ei kcatðiÞ :

ð10Þ

The implication of (9) for a two enzyme pathway is that if both enzymes approach saturation and changes in enzyme concentration are positive, then the sum of control coefficients approaches zero. In other words one can prove that if dEi > 0 for all i then limSat1 ;Sat2 !1 ðC1J þ C2J Þ ¼ 0:

ð11Þ

Proposition (11) holds due to the property that whenever two enzymes are near saturation, a situation is created whereupon the flux effects of increasing the concentration of one enzyme is eventually restrained by the other enzyme. This is because both enzymes are near their maximal rate and are hence rate limiting. Hence for any given positive change in enzyme concentration, the nearer the two enzymes approach saturation, the lower will be the sum of their control coefficients. When all changes in enzyme concentration are negative, one can obtain a reverse effect. For a two enzyme pathway 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. One can show that for an interval ai < Ei < bi ,   oSati ð12Þ 6 0 and ai P 0 and dEi < 0 for all i then limSat1 ;Sat2 !1 ðC1J þ C2J Þ ¼ 2: if Ei Proposition (12) illustrates a situation in which the pathway is not robust to decreases in enzyme concentration. The reason why (12) holds can be intuitively explained. For any given negative change in enzyme concentration, the control coefficient approaches one as the enzyme approaches

34

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

saturation and hence becomes rate limiting. If both enzymes approach saturation then their sum can exceed one. Eqs. (11) and (12) imply that for finite changes of any magnitude, the sum of control coefficients in a two enzyme pathway can range anywhere between 0 and 2. The implication is that dominance is not an inherent property of two enzyme pathways.

4. Sequential pathways of any length The two enzyme case serves as an initial guide to the properties of sequential pathways of length n > 2. These properties are a natural extension of the two enzyme case and we can summarize them in a brief format (see Appendix A for proofs). Consider the class of sequential two enzyme pathways that transform an input substrate s1 to a final product snþ1 and output o: Enzyme 1

Enzyme n

sink

s1 !    !    ! snþ1 ! o: Furthermore consider the case where the steady-state flux is given by a function gn : Rn 7! R1 such that ð13Þ J ¼ gn ðEÞ; where E ¼ hE1 ; . . . ; En i is the vector of total concentrations for enzymes 1 to n. 4.1. Epistasis Epistasis for a function of n variables can be measured by the off-diagonal entries of its Hessian matrix [74]. Let us define H as the Hessian for gn such that 1 0 2 oJ o2 J  B oE2 oE1 oEn C C B 1 C B C B .. .. .. ð14Þ H ¼B C: . . . C B C B 2 @ oJ o2 J A  oEn oE1 oEn2 One can show that the summation theorem implies that all off-diagonal entries of H are zero. Defining hjk as an individual entry in H , n X CiJ ¼ 1 for all dEi then hjk ¼ 0 for all j; k 2 f1; . . . ; ng where j 6¼ k: ð15Þ if i¼1

Proposition (15) P immediately implies that if any off-diagonal entry hjk 6¼ 0 then there exists a set of dEi where ni¼1 CiJ 6¼ 1. In addition one can show that all diagonal entries are also zero. In other words if

n X

CiJ ¼ 1 for all dEi then hjk ¼ 0 for all j; k 2 f1; . . . ; ng where j ¼ k:

i¼1

Proposition (16) is easier to prove than (15). In fact (16) is used to prove (15).

ð16Þ

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

35

4.2. Linearity Propositions (16) and (15) imply that the finite summation theorem holds if and only if gn is a linear function. In other words n n X X CiJ ¼ 1 if and only if J ¼ c i Ei ; ð17Þ i¼1

i¼1

where ci is a constant for every i. 4.3. Saturation With saturation one can prove that if dEi > 0 for all i 2 f1; . . . ; ng then for any ðj; lÞ 2 f1; . . . ; ng :

lim

Satj ;Satl !1

n X

! CiJ

¼ 0:

i¼1

ð18Þ Furthermore one can show that in the limiting case where all enzymes are approaching saturation for an interval ai < Ei < bi ,   oSati P 0 and ai P 0 and ! dEi < 0 for all i 2 f1; . . . ; ng if Ei n X ð19Þ then limSat1 ;...;Satn !1 CiJ ¼ n: i¼1

Eqs. (18) and (19) imply that for finite changes of any magnitude, the sum of control coefficients in a pathway with n enzymes can range anywhere between 0 and n.

5. Numerical analysis of a two enzyme pathway To illustrate the analytical results derived in the previous sections we examine a simple two enzyme pathway. A general dynamical model of a multi-enzyme system 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. The formulation of this pathway is inspired by the first three steps of the metabolism of lactose, which includes diffusion of lactose into the periplasmic space, active transport by b-galactoside permease and hydrolysis by b-galactosidase. This pathway has been extensively studied in Escherichia coli as a model system for metabolic evolution [59,60]. The scenario is represented by the following kinetic model shown schematically:

36

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51 D

u! D s1 ; k1

k3

k5

k7

! e1 þ s1 ! k2 ðes1 Þ k4 s2 þ e1 ;

ð20Þ

! e2 þ s2 ! k6 ðes2 Þ k8 s3 þ e2 ; Q

s3 !o; where ðes1 Þ and ðes2 Þ 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 (20) is a class of model that can be represented as a system of autonomous differential equations with input uðtÞ and initial conditions sð0Þ and eð0Þ such that s_ ðtÞ ¼ f1 ðsðtÞ; eðtÞ; uðtÞÞ; e_ ðtÞ ¼ f2 ðsðtÞ; eðtÞÞ; ð21Þ o_ ðtÞ ¼ f3 ðs3 ðtÞÞ; where at time t, the vector sðtÞ ¼ hs1 ðtÞ; s2 ðtÞ; s3 ðtÞi is the state vector of substrate concentrations, the vector eðtÞ ¼ he1 ðtÞ; es1 ðtÞ; e2 ðtÞ; es2 ðtÞi is the state vector of enzyme and enzyme–substrate complex concentrations, the variable oðtÞ is the concentration associated with the output variable and o_ ðtÞ is the output flux. The notation x_ associated with any variable x denotes the time derivative dx=dt. If the system reaches a steady-state flux the output flux o_ ðtÞ is equivalent to the steady-state flux J used in the MCA literature. During the interval of analysis 0 6 t 6 tf , the input remains constant such that 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 (21) combined with the kinetic model (20) are d½s1  ¼ D½s1   k1½e1 ½s1  þ D½u þ k2½es1 ; dt d½es1  ¼ k1½e1 ½s1   k2½es1   k3½es1  þ k4½e1 ½s2 ; dt d½s2  ¼ k3½es1   k4½e1 ½s2   k5½e2 ½s2  þ k6½es2 ; dt d½es2  ¼ k5½e2 ½s2   k6½es2   k7½es2  þ k8½e2 ½s3 ; dt d½s3  ¼ k7½es2   Q½s3   k8½e2 ½s3 ; dt d½e1  ¼ k1½e1 ½s1  þ k2½es1  þ k3½es1   k4½e1 ½s2 ; dt d½e2  ¼ k5½e2 ½s2  þ k6½es2  þ k7½es2   k8½e2 ½s3 ; dt d½o ¼ Q½s3 : dt

ð22Þ

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51 10

8

37

[E1] 6

4 2 0 1 0.75 0.5 J 0.25 0 10 8 6

4 2

[E2]

0

Fig. 1. Metabolic flux as a function of total enzyme concentrations. The constant environmental input is u ¼ 0:75 mM. The catalytic turnover rates are symmetric such that kcatð1Þ ¼ kcatð2Þ ¼ 42 000 min1 . Enzyme concentrations ½Ei  in lM. Flux in mM s1 .

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 where the derivatives s_ ðtÞ and e_ ðtÞ are equal to zero. Consider a case where the equations in (22) serve as the basis for the functions f1 , f2 and f3 with the stoichiometric constraint that ½E1  ¼ ½e1  þ ½es1 ; ½E2  ¼ ½e2  þ ½es2 :

ð23Þ

In such a case, the steady-state solution can be obtained in an implicit form (see Appendix B). The implicit equations can then be solved numerically by Newton–Raphson based algorithms that converge to the flux surface shown in Fig. 1. 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. 5.1. Test of the summation theorem We can test the summation of control coefficients P for the enzymatic steps directly. Figs. 2 and 3 show surfaces for for finite versions of the sum ni¼1 CiJ . It is evident from these examples that there is no evidence for an invariant relation in which the control coefficients sum to the predetermined value of one. Figs. 2 and 3 serve as illustrations of propositions (11) and (12). The regions where both enzymes approach saturation are shown in Fig. 4. Note that propositions (11) and (12) are independent of whether the pathway includes diffusion steps or not. The summation theorem fails in the regimes where both enzymes are near saturation or where flux is near its diffusion-limited maximal rate. The regions where the sum of control coefficients is equal to one are the regimes where flux is behaving linearly with respect to enzyme concentration.

38

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51 10 8 6

[E1] 4 2 0

1

0.8 0.6

J J C +C 1

2

0.4 0.2 0 10 8 6 4 2

[E2]

0

Fig. 2. Sum of control coefficients as a function of total enzyme concentrations. Enzyme concentrations ½Ei  in lM. For both enzymes dEi ¼ 0:2 lM. 5 4

[E1]

3 2

1 0 2 1.5 J

1

J

C 1+C 2

0.5 0

0

1

3

2

4

5

[E2]

Fig. 3. Sum of control coefficients as a function of total enzyme concentrations. Enzyme concentrations ½Ei  in lM. For both enzymes dEi ¼ 0:2 lM.

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 n X J J CiJ ¼ CDJ þ CE1 þ CE1 þ CQJ : ð24Þ i¼1

However this conjecture does not hold. Such a conjecture does eliminate the effects of the diffusion steps but the effect of the saturation steps will remain. Fig. 5 shows the failure of the modified summation conjecture.

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

39

A 1 0.8 0.6

Saturation E1

0.4 0.2 0 2

2 4

4 6

6

[E1]

8

[E2]

8 10

B 1 0.8 Saturation 0.6

E2

0.4 0.2 00

0 2

2 4

4 6

[E1]

6 8

[E2]

8 10

C 1 0.8 0.6

Saturation E1,E2

0.4 0.2 00

0 2

2 4

4 6

[E1]

6 8

8

[E2]

10

Fig. 4. Enzyme saturation as a function of total enzyme concentrations. Single enzyme saturation is measured as J =Ei kcatðiÞ . Double enzyme saturation is measured as ðJ =E1 k3 þ J =E2 k7 Þ=2. Enzyme concentrations ½Ei  in lM (keep in mind graphical repositioning of the origin when comparing to Fig. 2).

5.2. Epistatic interactions between enzymes The question of the extent to which enzymes have interactive effects on flux can also be addressed directly. In Fig. 6 the mixed partial difference d2 J =dE1 dE2 scaled to J is used as a measure of epistatic interactions between enzymes. Note that the region where the epistatic interactions increase coincides with a region where the summation theorem fails and both enzymes are approaching saturation.

40

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51 10 8 6

[E1] 4 2 0 1 0.8 0.6 C J+CJ+CJ+CJ D 1 2 Q 0.4 0.2 0 10

8 6 4 2

[E2]

0

Fig. 5. Failure of a modified summation conjecture. Enzyme concentrations ½Ei  in lM. For the control coefficients dEi ¼ 0:02 lM and dD=D ¼ dQ=Q ¼ 0:01. 10

8

[ E1 ] 6

4

2 0.05 0.04 0.03 1 2 . δ J 0.02 J δ Ε1 δΕ2 0.01 100 8

6 4 2

[ E2 ]

Fig. 6. Measure of epistatic interaction effects on flux. Enzyme concentrations ½Ei  in lM. For both enzymes dEi ¼ 0:02 lM.

6. Discussion and conclusion One of the central questions in evolutionary biology is how phenotypes evolve. The evolution of phenotypes in turn, depends on the evolution of development. Consequently to understand how phenotypes evolve we have to understand how development relates genotypes to phenotypes, and how this relationship can change throughout evolution. This is the context which compels us to examine the problem of phenotypic robustness. Can the robustness properties of a phenotype be altered throughout evolution. If so how? Admittedly this is a problem that goes beyond metabolic physiology. The mechanistic process that lead to robustness in different phenotypes do not solely lie within the realm of metabolic physiology. Along with the physical and chemical factors involved, phenotypes rely on an integration of metabolism, gene regulation and signal transduction. Hence metabolism is only a small piece of the puzzle. Nonetheless there are also generic properties

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

41

that have to be shared by all developmental processes that can undergo robustness modification. One of them is that robustness modification can occur as a result of mutations only if the genotype– phenotype relationship exhibits non-linearities in the form of epistasis. Hence when asking whether a given genotype–phenotype relationship can be modified for robustness, we should ask what kind of non-linear interactions the given relationship can exhibit. The intent of the present paper has been to present an analysis and critique of the biological conclusions derived from MCA with regard to phenotypic robustness in metabolism. Our results strongly support initial objections that were made with regards to the MCA perspective on the biochemical nature of dominance [29,51,65]. At the core of the problem of multi-enzyme systems is interaction effects between enzymes. In order to address these interactions, the non-linear properties of enzyme catalysis cannot be ignored. These interactions play a major role in the phenotypic characteristics of metabolic physiology and can lead to epistasis. Our analysis shows that for finite changes of any magnitude, the flux summation theorem does not hold for flux functions that exhibit any form of non-linearity. Since enzyme properties such as saturation do lead to nonlinearities and epistasis, it follows that there are no a priori constraints that would require the magnitude of mutational effects to be in a low range. This leaves open a possibility that had been previously rejected in MCA. Rate limiting steps or controlling enzymes do not have to be rare, and phenotypic robustness with respect to mutations is not an inevitable property of metabolic pathways. Furthermore propositions (18) and (19) suggest that the sum of control coefficients depends on the number of saturated enzymes. Since saturation is dependent on enzyme concentration and kinetic values such as kcat , it follows that robustness properties can be tuned through the modification of these values. This supports the position that there are no physiological impediments to the evolution of phenotypic robustness (or fragility for that matter) in metabolic pathways. As a guide for future work on the evolution of phenotypic robustness in metabolism there are two noteworthy caveats. First, our work does not consider the role that feedback regulation can play in the evolution of robustness [31]. Second, we do not consider the population dynamics under which robustness would evolve. The roles of regulation and population dynamics are matters that have to be addressed in future works on this topic. Acknowledgements Special thanks go to Michelle Girvan whose advice has been invaluable in improving the quality of this work. David KrakauerÕs advice greatly improved the final version of the manuscript. Ongoing discussions with Walter Fontana and Erica Jen on the topic of robustness are much appreciated. Leo Buss, Thomas Hansen, Michael Savageau and Lee Segel gave valuable advice on earlier incarnations of this work. Appendix A Proof of propositions (7) and (16). We start with the summation theorem definition of CiJ from (1) we know that n X dJi Ei ¼ 1 for all E 2 Rn and i 2 f1; . . . ; ng: dE J i i¼1

Pn

i¼1

CiJ ¼ 1. Using the

42

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

Hence J ¼ gn ðEÞ ¼

n X dJi Ei : dEi i¼1

ðA:1Þ

For any vector of length n let vi be a basis vector such that the ith term in vi is equal to unity and all other terms are zero. Hence n X Ei vi : ðA:2Þ E¼ i¼1

Now let us define a function dJi : Rnþ1 7! R1 such that for any E and dEi, dJi ðdEi ; EÞ  gn ðE þ dEi vi Þ  gn ðEÞ:

ðA:3Þ

For any arbitrary enzyme j 2 f1; . . . ; ng we can define fi 6¼ jg as the set of the n  1 enzymes other than j such that i 2 fi 6¼ jg if and only if i 2 f1; . . . ; ng and i 6¼ j. Let us choose dEj ¼ 1 for enzyme j and arbitrary dEi values for all other enzymes i 2 fi 6¼ jg such that 8i : dEi 6¼ 0. Using (A.1), for any given E we have X dJi ðdEi ; EÞ dJj ð1; EÞ J ¼ gn ðEÞ ¼ Ei : ðA:4Þ Ej þ 1 dEi fi6¼jg Similarly for the same E and the same set of dEi as in (A.4) we can choose any other dEj ¼ h such that h 6¼ 0 and X dJi ðdEi ; EÞ dJj ðh; EÞ J ¼ gn ðEÞ ¼ Ei : ðA:5Þ Ej þ h dEi fi6¼jg Subtracting (A.4) from (A.5) leads to the proposition that for every j and E, dJj ðh; EÞ ¼ hdJj ð1; EÞ:

ðA:6Þ

Separately consider two arbitrary positive real numbers h1 and h2 such that h ¼ h1 þ h2 . From the properties of real functions we know that dJj ðh; EÞ ¼ h1 dJj ð1; EÞ þ h2 dJj ð1; E þ h1 vj Þ:

ðA:7Þ

From equating (A.6) to (A.7) we have ðh  h1 ÞdJj ð1; EÞ ¼ h2 dJj ð1; E þ h1 vj Þ. But since h  h1 ¼ h2 , we are left with the proposition that for any dEj ¼ h1 , dJj ð1; EÞ ¼ dJj ð1; E þ dEj vj Þ:

ðA:8Þ

Note that by definition   o2 J dJj ðh; E þ hvj Þ  dJj ðh; EÞ ¼ lim : h2 oEj2 h!0 Hence (A.6) and (A.8) imply that for any E and j, o2 J ¼ 0: oEj2



ðA:9Þ

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

43

Proof of no epistasis when n ¼ 2 (proposition (5)). Consider a function gn that obeys the finite summation theorem (the function g being the particular case where n ¼ 2). Let E0 be the vector designating the origin such that Ei ¼ 0 for all i. Consider the case where n ¼ 2 and hi ¼ Ei . Let d Jei ðEÞ  dJi ð1; EÞ: ðA:10Þ By using (A.6) and the definitions in (A.1) and (A.3) we have J 1 ðE1 v1 þ E2 v2 ÞE1 þ d e J 2 ðE1 v1 þ E2 v2 ÞE2 : J ¼ gðE1 ; E2 Þ ¼ gn ðE1 v1 þ E2 v2 Þ ¼ d e

ðA:11Þ

Separately, by substituting Ei ¼ 0 for all i in (A.1) we know that gn ðE0 Þ ¼ 0:

ðA:12Þ

By (A.6) we also know that for any i and E, J i ðEÞ: dJi ðhi ; EÞ ¼ hi d e

ðA:13Þ

Hence by starting at the origin E0 , setting hi ¼ Ei for any i and making separate changes for E1 and E2 we can write J 1 ðE0 ÞE1 þ d e J 2 ðE1 v1 ÞE2 ðA:14Þ gn ðE1 v1 þ E2 v2 Þ ¼ d e and J 2 ðE0 ÞE2 þ d e J 1 ðE2 v1 ÞE1 : ðA:15Þ gn ðE1 v1 þ E2 v2 Þ ¼ d e J 1 ðE1 v1 þ E2 v2 Þ and d e J 2 ðE1 v1 Þ ¼ d e J 2 ðE1 v1 þ E2 v2 Þ. Hence But by (A.9) we know that d e J 1 ðE2 v2 Þ ¼ d e J 1 ðE0 ÞE1 þ d e J 2 ðE1 v1 þ E2 v1 ÞE2 ; gn ðE1 v1 þ E2 v2 Þ ¼ d e

ðA:16Þ

J 2 ðE0 ÞE2 þ d e J 1 ðE1 v1 þ E2 v1 ÞE1 : gn ðE1 v1 þ E2 v2 Þ ¼ d e

ðA:17Þ

Subtracting (A.16) from (A.11), for any E1 6¼ 0 and any E2 , we have de J 1 ðE1 v1 þ E2 v2 Þ ¼ d e J 1 ðE0 Þ:

ðA:18Þ

Similarly by substracting (A.17) from (A.11), for any E2 6¼ 0 and any E1 , we have de J 2 ðE1 v1 þ E2 v2 Þ ¼ d e J 2 ðE0 Þ:

ðA:19Þ

Note that by definition     o2 J dJ1 ðh; E þ hv2 Þ  dJ1 ðh; EÞ dJ2 ðh; E þ hv1 Þ  dJ2 ðh; EÞ ¼ lim ¼ lim : h!0 oE1 oE2 h!0 h2 h2

Hence (A.18) and (A.19) imply that for any E, o2 J ¼ 0: oE1 oE2



ðA:20Þ

Proof of no epistasis for any n > 2 (proposition (15)). For a pathway of length n > 2 consider a function P gn that obeys the finite summation theorem. We shall prove by induction that for any E ¼ ni¼1 Ei vi and any j 2 f1; . . . ; ng, ! n X J j ðE0 Þ: Ei vi ¼ d e ðA:21Þ de Jj i¼1

44

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

We first know that for any Ej 6¼ 0 and any Ek , where j, k 2 f1; . . . ; ng and j 6¼ k, J j ðE0 Þ: d Jej ðEj vj þ Ek vk Þ ¼ d e

ðA:22Þ The proof of (A.22) is virtually identical to the derivation of (A.18) with the only difference that n > 2 and (A.22) is valid for any enzymes j and k instead of just enzymes 1 and 2. Separately let us define a function Q that performs an arbitrary re-labelling of indices i 2 f1; . . . ; ng such that for every i there exists a unique ki 2 f1; . . . ; ng such that QðiÞ ¼ ki : ðA:23Þ Furthermore for any given Q let Q1 be the inverse function such that Q1 ðki Þ ¼ i. We can rewrite (A.22) in terms of an arbitrary choice of enzymes determined by Q such that for all Q, ! 2 X J Qð1Þ ðE0 Þ: EQðiÞ vQðiÞ ¼ d e ðA:24Þ de J Qð1Þ i¼1

For a proof by induction we need to prove for all n and all Q that ! m X if for all m; d JeQð1Þ J Qð1Þ ðE0 Þ where 2 6 m 6 p and p < n EQðiÞ vQðiÞ ¼ d e i¼1

then d e J Qð1Þ

pþ1 X

EQðiÞ vQðiÞ

! ¼ de J Qð1Þ ðE0 Þ:

ðA:25Þ

i¼1

Proof by induction. By the definition in (A.1) we know that for any E ¼ then we have ! ! n n n n X X X X dJi J ¼ gn Ei vi ¼ Ei ¼ de Ji Ei vi Ei : h i i¼1 i¼1 i¼1 i¼1 Hence for all Q, gn

n X

EQðiÞ vQðiÞ

i¼1

! ¼

n X

de J QðjÞ

j¼1

n X

Pn

i¼1

Ei vi if we let hi ¼ Ei

ðA:26Þ

! EQðiÞ vQðiÞ EQðjÞ :

ðA:27Þ

i¼1

In a similar fashion consider a case for which we pick p þ 1 enzymes for which hi ¼ Ei and Ei 6¼ 0 and n  p  1 enzymes for which Ei ¼ 0 and hi 6¼ 0. Let fQg be the set of all possible functions Q. Consider those Q 2 fQg for which any k 2 f1; . . . ; p þ 1g in the domain of Q is uniquely mapped to the index of one of the p þ 1 non-zero enzymes. It follows that ! ! pþ1 pþ1 pþ1 X X X ¼ E : gn E v d Je E v ðA:28Þ QðiÞ QðiÞ

i¼1

QðjÞ

j¼1

QðiÞ QðiÞ

QðjÞ

i¼1

Separately by using (A.12) and (A.13), if we start from the origin E0 , set hQðiÞ ¼ EQðiÞ and make separate changes for all EQðiÞ where i 2 f1; . . . ; p þ 1g we can write ! ! pþ1 pþ1 j1 X X X EQðiÞ vQðiÞ ¼ de J QðjÞ EQðiÞ vQðiÞ EQðjÞ : ðA:29Þ gn i¼1

j¼1

i¼1

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

45

P Note that in (A.29), when i ¼ 1 then i1 k¼1 dEQðkÞ vQðkÞ ¼ 0. But from (A.9) we also know that for any j, ! ! j1 j X X ðA:30Þ EQðiÞ vQðiÞ ¼ d JeQðjÞ EQðiÞ vQðiÞ : d JeQðjÞ i¼1

i¼1

Hence (A.29) becomes ! ! pþ1 pþ1 j X X X gn EQðiÞ vQðiÞ ¼ de J QðjÞ EQðiÞ vQðiÞ EQðjÞ : i¼1

j¼1

ðA:31Þ

i¼1

From the induction step in (A.25) we are assuming that for all m where 2 6 m 6 p and p < n, ! m X J QðiÞ ðE0 Þ: EQðiÞ vQðiÞ ¼ d e ðA:32Þ de J QðiÞ i¼1

The first p terms in (A.31) satisfy (A.32) so that we can rewrite (A.31) as ! ! pþ1 p pþ1 X X X EQðiÞ vQðiÞ ¼ de J QðjÞ ðE0 ÞEQðjÞ þ d e J Qðpþ1Þ EQðiÞ vQðiÞ EQðpþ1Þ : gn i¼1

j¼1

ðA:33Þ

i¼1

The last term in (A.33) is also the last term in (A.28), hence subtracting (A.33) from (A.28) we have ! p pþ1 p X X X E d Je E v ¼ d Je ðE0 ÞE : ðA:34Þ QðjÞ

j¼1

QðiÞ QðiÞ

QðjÞ

i¼1

QðjÞ

QðjÞ

j¼1

For the same set of p þ 1 enzymes we can choose p þ 1 different Q functions, each designated as Qk for each k 2 f1; . . . ; p þ 1g. Furthermore we can choose each Qk such that Qk ðp þ 1Þ ¼ k:

ðA:35Þ

Hence we can derive (A.34) for any k 2 f1; . . . ; pg and Qk . For each derivation using a given k and Qk , when equating the analogs of (A.33) and (A.28) the term corresponding to the d e J k function will be eliminated from the respective equation. Furthermore note that since all Q perform a one to one mapping from f1; . . . ; p þ 1g to the same set of p þ 1 enzymes, then for any two 0 Q; Q 2 fQg, pþ1 X

EQðiÞ vQðiÞ ¼

i¼1

pþ1 X

EQ0 ðiÞ vQ0 ðiÞ :

ðA:36Þ

i¼1

Thereby we have a system of p þ 1 linear equations which can be written in matrix form as Ax ¼ Ay where

ðA:37Þ

3  EQðiÞ vQðiÞ E1 7 6 7 6 .. x¼6 7 P .  5 4 pþ1 de J pþ1 E v E pþ1 i¼1 QðiÞ QðiÞ 2

de J1

P

pþ1 i¼1

ðA:38Þ

46

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

and 2

3 de J 1 ðE0 ÞE1 6 7 .. y¼4 5: . de J pþ1 ðE0 ÞEpþ1

ðA:39Þ

A is a ðp þ 1Þ  ðp þ 1Þ matrix for which Aij ¼ 1 if i 6¼ j and Aij ¼ 0 if i ¼ j. Since there is a unique zero term in each column of A, all the columns in A are linearly independent and the determinant satisfies det A 6¼ 0. Hence A1 exists, which means that x ¼ AA1 y and x ¼ y. Therefore for any k 2 f1; . . . ; p þ 1g, ! pþ1 X ¼ de J k ðE0 Þ: E v ðA:40Þ de Jk QðiÞ QðiÞ

i¼1

Eq. (A.40) can be derived for any set of p þ 1 enzymes. Hence (A.40) can be rewritten using any Q. Thereby proposition (A.25) is proven. Propositions (A.24) and (A.25) lead to (A.21) by induction. Note that by definition   o2 J dJj ðh; E þ hvi Þ  dJj ðh; EÞ ¼ lim : oEi oEj h!0 h2 Hence (A.21) implies that for any E and i; j 2 f1; . . . ; ng, o2 J ¼ 0: oEi oEj

ðA:41Þ



Example of induction when n ¼ 3. For the purpose of an illustration, here we use the induction proven for the general case n > 2 for the particular case where n ¼ 3. We know that Eq. (A.24) holds when p ¼ 2. Subsequently we only need one induction iteration since p þ 1 ¼ 3. Letting Q to perform the mappings 1 7! 1, 2 7! 2 and 3 7! 3, the analog of (A.28) is ! ! 3 3 3 X X X Ei vi ¼ d Jej Ei vi Ej : ðA:42Þ g3 i¼1

j¼1

i¼1

Let Q1 perform the mappings 1 7! 3, 2 7! 2 and 3 7! 1. As an analog of (A.29) we have ! 3 X J 3 ðE0 ÞE3 þ d e Ei vi ¼ d e J 2 ðE3 v3 ÞE2 þ d e J 1 ðE3 v3 þ E2 v2 ÞE1 : g3

ðA:43Þ

i¼1

Similarly letting Q2 do the mapping 1 7! 1, 2 7! 3 and 3 7! 2 we have ! 3 X J 1 ðE0 ÞE1 þ d e Ei vi ¼ d e J 3 ðE1 v1 ÞE3 þ d e J 2 ðE1 v1 þ E3 v3 ÞE2 : g3

ðA:44Þ

i¼1

Finally letting Q3 do the mapping 1 7! 1, 2 7! 2 and 3 7! 3 we have ! 3 X J 1 ðE0 ÞE1 þ d e Ei vi ¼ d e J 2 ðE1 v1 ÞE2 þ d e J 3 ðE1 v1 þ E2 v2 ÞE3 : g3 i¼1

ðA:45Þ

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

47

Using (A.9) and letting i 6¼ j 6¼ k we can substitute d Jei ðEj vj þ Ek vk Þ ¼ d e J i ðEi vi þ Ej vj þ Ek vk Þ in Eqs. (A.43), (A.44) and (A.45). We can subtract each of the resulting equations from (A.42), resulting in a system of three linear equations of the form Ax ¼ Ay as in (A.37), with 2 3 0 1 1 A ¼ 4 1 0 1 5: ðA:46Þ 1 1 0 The system of three linear equations leads to the solution that for all k 2 f1; 2; 3g and any E, ! 3 X Ei vi ¼ d Jek ðE0 Þ:  ðA:47Þ d Jek i¼1

Proof of linearity (propositions (8) and (17)). By definition   oJ hdJj ð1; EÞ ¼ lim : oEj h!0 h Propositions (A.9) and (A.41) imply that for any i there exists a constant ci such that for any E, oJ ¼ dJj ð1; EÞ ¼ ci : oEj

ðA:48Þ

From (A.1) we know that gn ðE0 Þ ¼ 0. Hence n X J ¼ gn ðEÞ ¼ ci Ei : i¼1

The reverse implication J ¼ tion. 

Pn

i¼1

ci Ei )

Pn

i¼1

ðA:49Þ CiJ ¼ 1 follows from direct algebraic manipula-

Proof of propositions (11) and (18). Consider an n enzyme pathway in which an enzyme j 2 f1; . . . ; ng is saturable such that for any E, J 6 Ej kcatðjÞ

and J 0 6 Ej kcatðjÞ ;

ðA:50Þ

where J 0 ¼ J þ dJ for any dEk where k 2 f1; . . . ; ng and k 6¼ j. It follows that if Ej is held constant, for any dEk > 0 where dJ P 0 we have Ej kcatðjÞ  J dJ 6 : dEk dEk

ðA:51Þ

Multiplying both sides by Ek =J we have Ek Ej kcatðjÞ  JEk Ek dJ 6 : J dEk J dEk

ðA:52Þ

Hence by (1) we obtain that when enzyme j is saturable, for any k 6¼ j, CkJ 6

Ek Ej kcatðjÞ Ek  : J dEk dEk

ðA:53Þ

48

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

Similarly if another enzyme l 6¼ j is saturable and El is held constant, for any m 6¼ l, we have Em El kcatðlÞ Em  : ðA:54Þ CmJ 6 J dEm dEm For each enzyme i, consider a measure of saturation Sati such that ðA:55Þ Sati ¼ J =Ei kcatðiÞ : It follows that for any enzyme i 2 fj; lg,   Ei kcatðiÞ lim ¼ 1: Sati !1 J Substituting (A.56) into (A.53) and (A.54) we have lim ðCkJ Þ ¼ 0

Satj !1

andSatl !1 ðCmJ Þ ¼ 0:

ðA:56Þ

ðA:57Þ

Given that j 6¼ l by definition, it follows that fk 2 f1; . . . ; ng j k 6¼ jg [ fm 2 f1; . . . ; ng j m 6¼ lg ¼ f1; . . . ; ng: Consequently using (A.57) and (A.58), when 8i dEi P 0 and j, l are saturable, then ! n X J  Ci ¼ 0: lim Satj ;Satl !1

ðA:58Þ

ðA:59Þ

i¼1

Proof of propositions (12) and (19). Consider an n enzyme pathway in which an enzyme j 2 f1; . . . ; ng is saturable such that for any dEj , J 6 Ej kcatðjÞ

and J 0 6 Ej0 kcatðjÞ ;

ðA:60Þ

where J 0 ¼ J þ dJ and Ej0 ¼ Ej þ dEj . Given that Sati ¼ J =Ei kcatðiÞ it follows that for any Ej , lim ðJ Þ ¼ Ej kcatðjÞ :

Satj !1

ðA:61Þ

Consider (A.60), (A.61) and the set of enzymes k such that k 2 f1; . . . ; ng and k 6¼ j. For domains in which all Ek are held constant and decreasing Ej increases saturation of enzyme j it follows that: For all intervals a < Ej < b, oSat2 if 6 0 and a P 0 and dEj < 0 then Ej lim ðdJ Þ ¼ Ej kcatðjÞ  Ej0 kcatðjÞ ¼ dEj kcatðjÞ :

Satj !1

Hence

 dEj kcatðjÞ =Ej kcatðjÞ dJ =J ¼ lim ¼ 1: ¼ lim Satj !1 Satj !1 dEj =Ej dE=E Similarly for any enzyme i that is saturable we can show that ðCjJ Þ

ðA:62Þ



lim ðCiJ Þ ¼ 1:

Sati !1

Hence for all intervals ai < Ei < bi ,

ðA:63Þ

ðA:64Þ

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

 if

oSati P 0 and ai P 0 and dEi < 0 Ei ! n X J lim Ci ¼ n:

Sat1 ;...;Satn !1

49

 for all i 2 f1; . . . ; ng then ðA:65Þ

i¼1

For a two enzyme pathway this means that   oSati if for all i; 6 0 ^ ai P 0 ^ dEi < 0 then Ei

lim

ðC1J þ C2J Þ ¼ 2:

Sat1 ;Sat2 !1



ðA:66Þ

Appendix B B.1. Derivations for numerical analysis 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: s1 ¼ ðe2 Qðk2 þ k3 Þk5 k7 þ e1 k2 k4 ðQðk6 þ k7 Þ þ e2 k6 k8 ÞÞH; es1 ¼ e1 k1 ðe2 Qk5 k7 þ e1 k4 ðQðk6 þ k7 Þ þ e2 k6 k8 ÞÞH; s2 ¼ e1 k1 k3 ðQðk6 þ k7 Þ þ e2 k6 k8 ÞH;

ðB:1Þ

es2 ¼ e1 e2 k1 k3 k5 ðQ þ e2 k8 ÞH; s3 ¼ e1 e2 k1 k3 k5 k7 H; where H¼

uD : e2 Qðk2 þ k3 Þk5 k7 D þ e1 ðQk2 k4 ðk6 þ k7 ÞD þ e2 ðQk1 k3 k5 k7 þ k2 k4 k6 k8 DÞÞ

ðB:2Þ

In addition recall the stoichiometric constraint that E1 ¼ e1 þ es1 ; E2 ¼ e2 þ es2 :

ðB:3Þ

Taking Eqs. (B.1)–(B.3) 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 fD; E1 ; E2 ; Q; k1 ; k2 ; k3 ; k4 ; k5 ; k6 ; k7 ; k8 g is given. The solution for fs1 ; s2 ; s3 ; e1 ; es1 ; e2 ; es2 g that is arrived at by numerical approximation can then be used to solve for the steady-state flux J . The surface in Fig. 1 is the numerical solution to (B.1)–(B.3). The wild-type kinetic parameters used in the numerical procedures were k1 ¼ k5 ¼ 4  107 1 1 M s , k2 ¼ k6 ¼ 4  102 s1 , k3 ¼ k7 ¼ 7  102 s1 and k4 ¼ k8 ¼ 1  106 M1 s1 . The diffusion constants used are D ¼ Q ¼ 3  101 s1 . For each enzyme Keq ¼ 70 and for the wild type kcat ¼ 42 000 min1 .

50

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

References [1] O. Mayo, R. Burger, Evolution of dominance: a theory whose time has passed?, Biol. Rev. 72 (1997) 97. [2] R.A. Fisher, The possible modification of the response of the wild type to recurrent mutations, Am. Natural. 62 (1928) 115. [3] W.L. Tower, The determination of dominance and the modification of behavior in alternative (Mendelian) inheritance, by conditions surrounding or incident upon the germ cells at fertilization, Biol. Bull. 18 (1910) 285. [4] C. Bridges, Non-disjunction of sex chromosome of Drosophila, J. Exp. Zool. 15 (1913) 587. [5] H.S. Jennings, Modifying factors and multiple allele morphs in relation to the results of selection, Am. Natural. 52 (1917) 301. [6] D.E. Lancefield, An autosomal bristle modifier affecting a sex-linked character, Am. Natural. 52 (1918) 462. [7] N.W. Timofeeff-Ressovsky, Studies on the phenotypic manifestation of hereditary factors. I. On the phenotypic manifestation of the genovariation radius incompletus in Drosophila funebris, Genetics 12 (1927) 128. [8] S. Wright, FisherÕs theory of dominance, Am. Natural. 63 (1929) 274. [9] J.B.S. Haldane, A note on FisherÕs theory of dominance, Am. Natural. 64 (1930) 87. [10] J.B.S. Haldane, The theory of the evolution of dominance, Am. Natural. 37 (1939) 365. [11] M.W. Feldman, S. Karlin, The evolution of dominance: a direct approach through the theory of linkage and selection, Theor. Popul. Biol. 2 (1971) 482. [12] C.A. Clarke, P.M. Sheppard, The evolution of dominance under disruptive selection, Heredity 14 (1960) 73. [13] P. OÕDonald, J.A. Barrett, Evolution of dominance in polymorphic Batesian mimmicry, Theor. Popul. Biol. 4 (1973) 173. [14] D. Charlesworth, B. Charlesworth, Theoretical genetics of Batesian mimicry III. The evolution of dominance, J. Theor. Biol. 55 (1975) 325. [15] R. Burger, Nonlinear analysis of some models for the evolution of dominance, J. Math. Biol. 16 (1983) 269. [16] S.P. Otto, D. Bourguet, Balanced polymorphism and the evolution of dominance, Am. Natural. 153 (1999) 561. [17] H. Kacser, J.A. Burns, The molecular basis of dominance, Genetics 97 (1981) 639. [18] R. Burger, On the evolution of dominance modifiers. I. A nonlinear analysis, J. Theor. Biol. 101 (1983) 585. [19] R. Burger, Dynamics of the classical genetic model for the evolution of dominance, Math. Biosci. 67 (1983) 125. [20] G.P. Wagner, R. Burger, A non-equilibrium approach to the evolution of genetic systems, J. Theor. Biol. 113 (1985) 475. [21] H. Kacser, J.A. Burns, D. Fell, The control of flux: 21 years on. The control of flux, Bioch. Soc. Trans. 23 (1995) 341. [22] H. Kacser, The control of flux: 21 years on. Recent developments in Metabolic Control Analysis, Bioch. Soc. Trans. 23 (1995) 387. [23] H.A. Orr, A test of FisherÕs theory of dominance, Proc. Nat. Acad. Sci. USA 88 (1991) 11413. [24] P.D. Keightley, A metabolic basis for dominance and recessivity, Genetics 143 (1996) 621. [25] J.W. Porteous, Dominance-one hundred and fifteen years after MendelÕs paper, J. Theor. Biol. 182 (1996) 223. [26] C.D. Meiklejohn, D.L. Hartl, A single mode of canalization, Trend. Ecol. Evol. 17 (2002) 468. [27] J. Hodgkin, Fluxes, doses and poisons: molecular perspectives on dominance, Trends Genet. 9 (1993) 1. [28] D.R. Forsdyke, The heat–shock response and the molecular basis of genetic dominance, J. Theor. Biol. 167 (1994) 1. [29] M.A. Savageau, Dominance according to metabolic control analysis: major achievement or house of cards?, J. Theor. Biol. 154 (1992) 131. [30] U. Grossniklaus, M. Madhusudhan, V. Nanjundiah, Nonlinear enzyme kinetics can lead to high metabolic flux control coefficients: implications for the evolution of dominance, J. Theor. Biol. 182 (1996) 299. [31] S. Omholt, E. Plahte, L. Oyehaug, K. Xiang, Gene regulatory networks generating the phenomena of additivity, dominance and epistasis, Genetics 155 (2000) 969. [32] M.A. Gilchrist, H.F. Nijhout, Nonlinear developmental processes as sources of dominance, Genetics 159 (2001) 423. [33] D. Bourguet, The evolution of dominance, Heredity 83 (1999) 1. [34] C.H. Waddington, The strategy of the genes, MacMillan, New York, NY, 1957. [35] J.M. Rendel, Canalization and gene control, Academic Press, New York, NY, 1967. [36] I.I. Schmalhausen, Factors of evolution: the theory of stabilizing selection, University of Chicago Press, Chicago, IL, 1986. [37] G.P. Wagner, G. Booth, H. Bagheri-Chaichian, A population genetic theory of canalization, Evolution 51 (1997) 329.

H. Bagheri-Chaichian et al. / Mathematical Biosciences 184 (2003) 27–51

51

[38] S.C. Stearns, M. Kaiser, T.J. Kawecki, The differential canalization of fitness components against environmental perturbations in Drosophila melanogaster, J. Evol. Biol. 8 (1995) 539. [39] A. Wagner, Does evolutionary plasticity evolve?, Evolution 50 (1996) 1008. [40] N. Barkai, S. Leibler, Robustness in simple biochemical networks, Nature 387 (1997) 913. [41] J. Gerhart, M. Kirschner, Cells, Embryos and Evolution, Blackwell Science, Malden, MA, 1997. [42] A. Wagner, Robustness against mutations in genetic networks of yeast, Nat. Genet. 24 (2000) 355. [43] G. Gibson, G.P. Wagner, Canalization in evolutionary theory: a stabilizing theory, Bioessays 22 (2000) 372. [44] L.W. Ancel, W. Fontana, Plasticity, evolvability and modularity in RNA, J. Exp. Zool. (Mol. Dev. Evol.) 288 (2000) 242. [45] J.L. Hartman, B. Gavrik, L. Hartwell, Principles for the buffering of genetic variation, Science 291 (2001) 1001. [46] D.C. Krakauer, J.B. Plotkin, Redundancy, antiredundancy, and the evolution of genomes, Proc. Nat. Acad. Sci. USA 99 (2002) 1405. [47] H. Kacser, J.A. Burns, The control of flux, Symp. Soc. Exp. Biol. 27 (1973) 65. [48] R. Heinrich, T.A. Rapoport, A linear steady-state treatment of enzymatic chains. General properties, control and effector strength, Eur. J. Biochem. 42 (1974) 89. [49] C. Reder, Metabolic control theory: a structural approach, J. Theor. Biol. 135 (1988) 175. [50] M.A. Savageau, Biochemical Systems Analysis: a study of function and design in molecular biology, AddisonWesley, Reading, MA, 1976. [51] M.A. Savageau, A. Sorribas, Constraints among molecular and systemic properties: implications for physiological genetics, J. Theor. Biol. 141 (1989) 93. [52] S. Schuster, R. Heinrich, The definitions of metabolic control analysis, Biosystems 27 (1992) 1. [53] R. Heinrich, S. Schuster, The regulation of cellular systems, Chapman and Hall, New York, NY, 1996. [54] P.D. Keightley, H. Kacser, Dominance, pleiotropy and metabolic structure, Genetics 117 (1987) 319. [55] P.D. Keightley, Metabolic models of selection response, J. Theor. Biol. 182 (1996) 311. [56] D. Fell, Understanding the Control of Metabolism, Portland, London, 1997. [57] D.L. Hartl, D. Dykhuizen, A.M. Dean, Limits of adaptation: the evolution of selective neutrality, Genetics 111 (1985) 655. [58] A.M. Dean, D. Dykhuizen, D.L. Hartl, Fitness as a function of b-galactosidase activity in Escherichia coli, Genet. Res. 48 (1986) 1. [59] D. Dykhuizen, A. Dean, D.L. Hartl, Metabolic flux and fitness, Genetics 115 (1987) 25. [60] A.M. Dean, Selection and neutrality in lactose operons of Escherichia coli, Genetics 123 (1989) 441. [61] D.L. Hartl, C.H. Taubes, Compensatory nearly neutral mutations: selection against adaptation, J. Theor. Biol. 182 (1996) 303. [62] R. Small, H. Kacser, Responses of metabolic systems to large changes in enzyme activities and effectors. 1. The linear treatment of unbranched chains, Eur. J. Biochem. 213 (1993) 613. [63] H. Kacser, H.M. Sauro, L. Acarenza, Enzyme–enzyme interactions and control analysis. 1. The case of additivity: monomer oligomer associations, Eur. J. Biochem. 187 (1990) 481. [64] B. Kholodenko, D. Molenaar, S. Schuster, R. Heinrich, H. Westerhoff, Defining control coefficients in non-ideal metabolic pathways, Biophys. Chem. 56 (1995) 215. [65] A. Cornish-Bowden, Dominance is not inevitable, J. Theor. Biol. 125 (1987) 333. [66] H. Kacser, Dominance not inevitable but very likely, J. Theor. Biol. 126 (1987) 505. [67] H. Kacser, A superior theory?, J. Theor. Biol. 149 (1991) 141. [68] B. Kholodenko, J. Rohwer, M. Cascante, H. Westerhoff, Subtleties in control by metabolic channelling and enzyme organization, Mol. Cell. Biochem. 184 (1998) 311. [69] T. Hofer, R. Heinrich, A second order approach to metabolic control analysis, J. Theor. Biol. 164 (1993) 85. [70] H. Kacser, L. Acarenza, A universal method for achieving increases in metabolite production, Eur. J. Biochem. 216 (1993) 361. [71] L. Acarenza, Design of large metabolic responses, constraints and sensitivity analysis, J. Theor. Biol. 207 (2000) 265. [72] C. Giersch, Control analysis of metabolic networks, Eur. J. Biochem. 174 (1988) 509. [73] S. Rice, The evolution of canalization and the breaking of Von BaerÕs law: modeling the evolution of development with epistasis, Evolution 52 (1998) 647. [74] S.H. Rice, The evolution of developmental interactions, in: J.B. Wolf, E.D. Brodie, M.J. Wade (Eds.), Epistasis and the Evolutionary Process, Oxford Unicversity, New York, NY, 2000, p. 82.