We have determined the effects of control by overall
feedback inhibition on the systemic behavior of unbranched metabolic pathways with an arbitrary pattern of other feedback inhibitions by
using a recently developed numerical generalization of Mathematically Controlled Comparisons, a method for comparing the function of alternative molecular designs. This method allows the rigorous determination of the changes in systemic properties that can be exclusively attributed to overall feedback inhibition. Analytical results show that the unbranched pathway can achieve the same steady-state flux, concentrations, and logarithmic gains with respect
to changes in substrate, with or without overall feedback inhibition.
The analytical approach also shows that control by overall feedback
inhibition amplifies the regulation of flux by the demand for end
product while attenuating the sensitivity of the concentrations to the
same demand. This approach does not provide a clear answer regarding
the effect of overall feedback inhibition on the robustness, stability,
and transient time of the pathway. However, the generalized numerical
method we have used does clarify the answers to these questions. On
average, an unbranched pathway with control by overall feedback
inhibition is less sensitive to perturbations in the values of the
parameters that define the system. The difference in robustness can
range from a few percent to fifty percent or more, depending on the length of the pathway and on the metabolite one considers. On average,
overall feedback inhibition decreases the stability margins by a
minimal amount (typically less than 5%). Finally, and again on
average, stable systems with overall feedback inhibition respond faster
to fluctuations in the metabolite concentrations. Taken together, these
results show that control by overall feedback inhibition confers
several functional advantages upon unbranched pathways. These
advantages provide a rationale for the prevalence of this control
mechanism in unbranched metabolic pathways in vivo.
 |
INTRODUCTION |
Biochemical control systems have been studied for
more than 45 years. The discovery of control by molecular feedback
inhibition in biochemical pathways was initially made in unbranched
biosynthetic pathways (Umbarger, 1956
; Yates and
Pardee, 1956). In these pathways, the most common pattern of
control is inhibition of the initial reaction by the final product of
the pathway (end-product inhibition or overall feedback inhibition).
There are several criteria for the functional effectiveness of control
in such pathways that can be used to evaluate the biological significance of the overall feedback inhibition mechanism. A
biochemical pathway should be robust, i.e., it should function
reproducibly despite perturbations in the values of the parameters that
define the structure of the system. The operating point (state) of the system should be stable so that the system returns to the steady state
following small random fluctuations in the values of the dependent
variables; if not, the system tends to be dysfunctional because
spurious environmental fluctuations will lead to loss of the steady
state. The flux through the pathway should be responsive to changes in
the demand for the final product. This ensures that the amount of
material flowing through the pathway is intimately coupled to the
metabolic needs of the cell. Finally, the system should be temporally
responsive to changes, because, otherwise, the system is unlikely to be
competitive in rapidly changing environments. [A more extensive
discussion of these criteria and their quantification can be found in
Savageau (1976) and Hlavacek and Savageau
(1997).]
There have been several studies focused on the effect of control by
overall feedback inhibition on the stability of unbranched pathways. In
general, the first enzyme of the pathway is considered to be
allosteric, whereas the others are considered to be Michaelian (e.g.,
Goodwin, 1963; Morales and McKay, 1967;
Walter, 1969a,b, 1970; Viniegra-Gonzalez,
1973; Hunding, 1974; Rapp, 1976;
Dibrov et al., 1981). The stability of an unbranched
pathway with overall feedback inhibition and enzymes confined to one of
two spatial compartments with diffusion between compartments has been
studied by Costalat and Burger (1996). They found that
stability can be increased by this type of compartmentation. These
studies considered pathways with no internal feedback inhibitions.
Several other patterns involving control by inhibitory feedback can, in
principle, perform the same qualitative functions as overall feedback
inhibition. One such pattern is, for example, a sequence of feedback
inhibitions in which each intermediate inhibits the reaction that
immediately precedes it (Koch, 1967). Other patterns of
internal feedback inhibition can be found by searching either the
literature or some of the databases for metabolism that are burgeoning
on the world wide web (e.g., KEEG:
http://www.genome.ad.jp/kegg/; ECOCYC:
http://ecocyc.PangeaSystems.com/ecocyc/server.html; PUMA: http://www.unix.mcs.anl.gov/compbio/PUMA/Production/puma_graphics.html; EMP: http://wit.mcs.anl.gov//EMP/). However, even when
intermediate feedback inhibition patterns exist, control by overall
feedback inhibition remains a prevalent theme in biosynthetic pathways.
Savageau (1972, 1974,
1975, 1976) studied the function of various patterns of feedback
inhibition and explained the prevalence of control by overall feedback
inhibition by using arguments based on selection. He assumed that the
design of a pathway is selected to optimize certain systemic
characteristics, and then compared those systemic characteristics in
unbranched pathways with overall feedback inhibition to the same
characteristics in pathways with alternative inhibitory feedback
designs. He showed that the pathway with control by overall feedback
inhibition is more robust, i.e., less sensitive to perturbations in
parameter values than the pathway with many alternative designs
(Savageau, 1974).
The stability of cases with control by internal feedback inhibitions
has also been examined (e.g., Savageau, 1976;
Thron, 1991a,b;
Demin and Kholodenko, 1993). These authors found that systems with internal feedback inhibitions have larger stability margins than systems without these interactions. They also determined that, for systems without internal feedback inhibition, control by
overall feedback inhibition decreases the stability margins of the pathway.
In this paper, we consider unbranched pathways with all possible
patterns of internal feedback inhibitions (the "fully-wired" case)
and use all of the criteria mentioned above to determine the biological
significance of control by overall feedback inhibition in such
pathways. We use a technique called Mathematically Controlled Comparison that was originally developed to determine irreducible qualitative differences in systemic behavior of models with alternative regulatory designs for the same network of reactions (Savageau, 1972, 1976; Irvine
and Savageau, 1985). This qualitative technique requires the
existence of closed-form solutions for the steady state. Such solutions
can be obtained by using the local S-system representation to
characterize the pathway of interest. Important functional constraints
are introduced by equating relevant steady-state properties of the
alternative systems being compared. The limitations of this technique
have been overcome by a recently developed generalization that uses
numerical methods to obtain results that are general in a statistical
sense (Alves and Savageau, 2000a).
| |
METHODS |
Alternative models and key systemic properties
Consider the unbranched pathways depicted in Fig.
1. The independent variable
Xn+1 represents the cell's demand for the end
product Xn. If the cell requires large amounts
of Xn, then the value of
Xn+1 will be high; if small amounts of
Xn are required, then the value of
Xn+1 will be low. The dynamic behavior of such
systems can be described in principle by a set of ordinary differential
equations. There is no generic representation of these equations that
can provide a globally accurate description of the behavior [see
Appendix]. However, the set of equations can be approximated to the
first order in logarithmic space (Savageau, 1969),
yielding ordinary differential equations with the canonical form of an
S-system (Savageau, 1996). This representation has a
solid theoretical foundation based on Taylor's theorem. Thus, the
validity of the results is guaranteed within some neighborhood of the
nominal steady-state operating point. The size of this neighborhood
cannot be specified in general, because it depends on the
characteristics of each individual system.
View larger version (12K):
[in this window]
[in a new window]
|
FIGURE 1
(A) Model of an unbranched pathway with all
possible inhibitory feedback interactions (reference model).
(B) Model of an unbranched pathway with all possible
inhibitory feedback interactions except overall feedback inhibition
(alternative model). The horizontal arrows represent biochemical
reactions, whereas the vertical arrows represent inhibitory feedback
interactions.
|
|
For pathways with n intermediates, the general case in which
all possible feedback inhibitions exist (Fig. 1 A) can be
described in the local S-system representation as
|
(1)
|
|
(2)
|
|
(3)
|
The corresponding case without overall feedback inhibition (Fig.
1 B) can be described by the same set of equations, except that Eq. 1 is replaced by
|
(4)
|
The rate law for each reaction is represented by a simple
product of power-law functions. The values for the parameters in this
representation can be determined directly from conventional experimental measurements of initial rate as a function of reactant and
modifier concentrations (Savageau, 1976). The range of
values for the concentrations is chosen to sample the region about the nominal steady state of interest.
The parameters are defined according to Taylor's theorem as
|
(5)
|
|
(6)
|
where the additional subscript zero signifies that the variables
and their derivatives are evaluated at the steady-state operating
point. The definition of these parameters allows them to be directly
related to the parameters in other representations such as the
traditional Michaelis-Menten representation. In the simplest case of
the Hill rate law,
|
(7)
|
[and the irreversible Michaelis-Menten rate law
(n = 1)], these relationships are well known
(Savageau, 1971a),
|
(8)
|
|
(9)
|
The multiplicative parameters, 
, can be interpreted as rate
constants that are always positive. The exponential parameters, g, can be interpreted as kinetic orders that represent the
direct influence of each intermediate on each rate law. If
Xi is directly involved in the rate law
Vj, either as a substrate or a modulator, and if
an increase in Xi causes an increase in the rate
Vj, then the kinetic order will be positive. If
an increase in Xi causes a decrease in
Vj, then the kinetic order will be negative. If Xi is not directly involved in
Vj, then the kinetic order will be zero. The
positive kinetic orders in Eqs. 1-4 are
gi+1,i (0 
i n) and
g'10, because these are the kinetic
orders for substrates of reactions, and
gn+1,n+1, which, together with
Xn+1, represents the demand for the end product
Xn. The remaining kinetic orders, which
represent feedback inhibitions, are negative.
At a steady state, the rate of production and the rate of consumption
will be equal for each intermediate, and Eqs. 1-3 reduce to the
following matrix equation (Savageau, 1969), which can be solved analytically.
|
(10)
|
where b1 = log(2/1),
bi = log(i+1/i),
aij = gij 
gi+1,j, and Yi = log(Xi). Eq. 10 is linear and therefore easily solved to obtain the steady-state value for each
Yi, and then the corresponding value for each
Xi is obtained by simple exponentiation. Eqs.
2-4 reduce to an identical matrix equation, except that the parameters
of the first row are primed and g'1n = 0.
Two types of systemic coefficients, logarithmic gains and parameter
sensitivities, can be used to characterize the steady state of such
models (Shiraishi and Savageau, 1992). Logarithmic gains
measure the relative influence of each independent variable on each
dependent variable of the integrated model. For example,
|
(11)
|
measures the percent change in the concentration of intermediate
Xi caused by a percentage change in the
concentration of the initial substrate X0.
Logarithmic gains provide important information concerning the
amplification or attenuation of signals as they are propagated through
the system. The experimental measurement of a logarithmic gain involves
the determination of steady-states fluxes and concentrations at
different values for a given independent variable (Savageau,
1971a).
Parameter sensitivities measure the relative influence of each
parameter on each dependent variable of the model. For example,
|
(12)
|
measures the percent change in the concentration of intermediate
Xi caused by a percentage change in the value of
the parameter pj. Parameter sensitivities
provide important information about system robustness, i.e., how
sensitive the system is to perturbations in the parameters that define
the structure of the system. Because enzymes usually have a first-order
influence on the process they catalyze, the logarithmic gain in flux
and in each concentration with respect to change in the concentration
of each enzyme is the same as the sensitivity in flux and in each
concentration with respect to change in the rate constant of the
corresponding enzyme. The experimental measurement of a parameter
sensitivity involves the determination of steady-state fluxes and
concentrations before and after changing the value of a parameter by
mutation or other means (Savageau, 1971b).
Because we can calculate closed-form steady-state solutions for Eqs.
1-3 and 2-4, we can also calculate each of the two types of
coefficients simply by taking the appropriate derivatives of those
solutions. Although the mathematical operations involved are the same
in each case, it is important to keep in mind that the biological
significance of the two types of coefficients is very different.
The local stability of the steady state can be determined by applying
the Routh criteria (Dorf, 1992). The magnitude of the two critical Routh conditions can be used to quantify the margin of
stability (e.g., Savageau, 1976).
The use of the S-System formalism allows for an analytical study of the
dynamical systems at steady state. Comparisons of systems with only one
feedback inhibition to systems without feedback regulation can be done
and interpreted in a fully symbolic way. However, for comparisons
involving many feedback inhibitions, numerical values must be
introduced for the parameters to make the comparisons interpretable.
The steady-state behavior of the alternative models is compared with
respect to their flux, intermediate concentrations, logarithmic gains
with respect to changes in initial substrate and demand for end
product, robustness, and stability margins. The differential equations
also are solved numerically to characterize the temporal responsiveness
of the alternative designs. To evaluate this, we increase the
steady-state concentration of each Xi by 20%
and measure the time the system takes to relax back to within 1% of
its original steady state.
Calculating constraints for the mathematically controlled
comparison
Only the first step in the pathway is allowed to differ between
the reference model (Fig. 1 A) and the alternative model
(Fig. 1 B). Therefore, to establish "internal
equivalence" (Savageau, 1972, 1976; Irvine, 1991) between the two designs,
we require the values for the corresponding parameters of all other
steps in the two models to be the same.
The first step of the pathway differs between the reference model and
the alternative model, and the degrees of freedom associated with this
difference must be eliminated to the extent possible. If we reason that
loss or gain of an inhibitory site on the first enzyme comes about by
mutation, and that this mutation can cause changes in all the
parameters of the first reaction, then a mutation causing loss of
overall feedback inhibition would change the parameters 1 and g10 through
g1n in Eq. 1 to the corresponding primed
parameters in Eq. 4. Clearly, the value of the parameter
g'1n, which equals zero, differs from
that of g1n, which is nonzero. The remaining primed parameters also will have values that, in general, are not equal
to the values of the corresponding parameters in the reference model.
Because we wish to determine those effects that are due solely to
changes in the structure of the system and not simply to arbitrary
changes in the values of parameters, we shall specify values for the
primed parameters that minimize all other effects. This can be
accomplished by deriving the mathematical expression for a given
steady-state property in each of the two models, equating these
expressions, and then solving the constraint equation for the value of
a primed parameter. This process establishes an "external
equivalence" between the two designs (Savageau, 1972, 1976; Irvine,
1991). After values for all the primed parameters have been
specified in terms of the known values for the reference system, the
extra degrees of freedom have been eliminated, and we can proceed with
the comparison.
Three classes of constraint equations are used to fix the values for
the k + 2 primed parameters when there are k
interactions that feed back to the first step of the alternative
pathway. These are obtained by equating steady-state logarithmic gains,
concentrations, and parameter sensitivities as described below.
First, equating the logarithmic gains for any one of the metabolites
with respect to change in the initial substrate,
|
(13)
|
which causes each of the other corresponding intermediates to
have the same logarithmic gain, specifies the value of the kinetic
order g'10 in terms of known values for
the reference system. This condition also makes the corresponding
logarithmic gain in flux the same for the two designs.
Second, equating the concentrations for any one of the metabolites in
the pathway,
|
(14)
|
which causes each of the corresponding intermediates to have the
same concentration, specifies the value of the rate constant '1. This condition also makes the flux the same
for the two designs.
Finally, the remaining k 1 primed parameters are
fixed by equating the rate-constant parameter sensitivities,
|
(15)
|
for any Xi and k 1
different rate constants j. Different results will be
obtained, depending upon which of the parameter sensitivities are not
used in this procedure.
For example, consider the case in which all n 1
intermediates feed back on the first step in the pathway. If the
unconstrained sensitivity in Eq. 15 is
S(Xn, n), then the values of the
primed parameters are given by
|
(16)
|
|
(17)
|
|
(18)
|
If the unconstrained sensitivity in Eq. 15 is
S(Xi, 1), then the values of the
primed parameters are
|
(19)
|
|
(20)
|
If the unconstrained sensitivity in Eq. 15 is
S(Xi, j) where 1 < j < n, then the values of the primed parameters are
|
(21)
|
|
(22)
|
|
(23)
|
Because the objective of a controlled comparison is to minimize
the differences between the systems being compared, we chose the
unconstrained sensitivity that leads to the smallest number of systemic
properties with values that differ between the reference system and the
alternative system. The systemic differences are minimized when the
unconstrained sensitivity is
S(Xi, n+1); any other choice
leads to at least one additional systemic property that differs between
the two systems.
If only a subset of the intermediates feed back on the first step of
the pathway, and if we use the constraint set that causes the smallest
number of properties to be different between systems A and B, then each
kinetic order representing a feedback inhibition has the same value in
both models, except for the kinetic order representing the last
intermediate to feed back on the first step of the pathway. In general,
|
(24)
|
where Xk is the last intermediate to feed
back on the first step of the pathway, and 
nk is a
positive subdeterminant of [A] that depends on the actual
Xk and on the length n of the
pathway. The kinetic orders g1p with
p < k are the same for both systems. As for the rate
constant '1, its general form is
|
(25)
|
where 
kp is either a function of the kinetic
orders or zero.
For the special case in which the final product is the only metabolite
to feed back on the initial step, the primed parameters are given by
|
(26)
|
|
(27)
|
This means that g'10 is always
smaller than g10. (To contrast these results
with the analogous results expressed within the Michaelis-Menten
formalism, see the Appendix.)
Numerical comparison
It is straightforward to compare analytically corresponding
magnitudes from each of the two designs. For two- and three-step pathways, the comparisons are clearly interpretable for most systemic properties. The analytical results give qualitative information that
characterizes the role of overall feedback inhibition for the system in
Fig. 1 A. As the length of the pathway increases, the
analytical interpretation becomes problematic. To determine if a given
magnitude is larger in the reference system or the alternative system
requires knowledge of actual parameter values in these cases and
a method, such as that found in Alves and Savageau (2000a), for making the numerical equivalent of a
mathematically controlled comparison.
To obtain numerical information, one must introduce specific values for
the parameters and compare systems. For this purpose, we have randomly
generated a large ensemble of parameter sets and selected 5000 of these
sets that define systems consistent with various physical and
biochemical constraints. These constraints include mass balance, low
concentrations of intermediates and small changes in their values to
minimize utilization of the limited solvent capacity in the cell, small
values for parameter sensitivities so as to desensitize the system to
spurious fluctuations affecting its structure, and stability margins
large enough to ensure local stability of the systems. A detailed
description of these methods can be found in Alves and Savageau
(2000c). Mathematica (Wolfram, 1997) was used
for all the numerical procedures. Pathways of up to seven steps were
studied using this numerical methodology.
To interpret the ratios that result from our analysis, we use density
of ratios plots as defined in Alves and Savageau
(2000b). The primary density plots of the raw data have the
magnitude of some property for the reference system on the
x-axis and the corresponding ratio of magnitudes
(alternative system to reference system) on the y-axis. The
primary plot can be viewed as a list of 5000 paired values that can be
ordered with respect to the reference magnitude, thereby forming a list
L1 in which the first pair has the lowest measured value for property P in the reference model, the
second has the second lowest, and so on.
Secondary density plots are constructed from the primary plots by the
use of moving quantile techniques with a window size of 500. The
procedure is as follows. One collects the first 500 ratios from the
list L1, calculates the quantile of interest for this sample, and pairs this number 
R
with the median
value of the corresponding P values for the reference model,
denoted P. One advances the window by one position,
collects ratios 2-501, calculates R, and pairs it
with the corresponding P value and continues in this
manner until the last ratio from the list L1 was
used for the first time (for further explanation of moving median
techniques see, e.g., Hamilton, 1994).
The slope in the secondary plot measures the degree of correlation
between the quantities plotted on the x- and
y-axes. This technique also is used to examine correlations
between ratios of interest and other magnitudes shared by the two
systems, e.g., the correlation between the ratio of stability margins
and the magnitude of a rate constant common to the two systems (for
traditional applications of correlation analysis, see Wherry,
1984).
| |
RESULTS |
Mathematically controlled comparison
Response to availability of substrate and demand for end
product
The responsiveness of each system to changes in the independent
concentration variables X0, which represents the
availability of initial substrate, and Xn+1,
which represents the demand for end product, is characterized by a set
of logarithmic gains that provides a quantitative measure of signal
propagation through the system.
The logarithmic gains of the two systems in response to changes in the
initial substrate are identical at each step in the pathway [i.e.,
L(Vi, X0)A = L(Vi, X0)B and
L(Xi, X0)A = L(Xi, X0)B for
1 i n] because of the constraints for
external equivalence described in the Methods section. Hence, the
responsiveness of the two systems to changes in the availability of
initial substrate is identical.
In contrast, the responsiveness of the two systems to changes in the
demand for their end product is different. The ratio of the logarithmic
gains in flux is given by
|
(28)
|
where 
is always a negative sum of products of the kinetic
orders, g1n < 0, and
gj+1,j > 0 for j = 1, 2, ... , n 1. These results demonstrate that the flux in
the reference system is more responsive than that in the alternative
system to changes in demand for end product.
The ratio of the logarithmic gains in concentration is given by
|
(29)
|
where is a sum of products of the kinetic orders that depends
on i and the length of the pathway,
g1n < 0, and
gj+1,j > 0 for j = 1, 2, ... , n 1. When i = 1 or i = n, is always positive and, thus, the reference model is
always less sensitive to demand. When 1 < i < n,
is positive in most cases. This shows that the concentrations are
usually less sensitive to demand in the system with overall feedback inhibition.
Robustness of flux
The robustness of any systemic property with respect to
perturbations in the values of the parameters that define the system is
characterized by a set of parameter sensitivities. The steady-state flux of reference and alternative systems has different sensitivities with respect to the parameters n, n+1,
g1,n
1, gn+1,n, gnn, and gn,n1 that are
common to the two systems. The sensitivities are the same with respect
to all other parameters common to the two systems.
The sensitivities of the steady-state flux with respect to the
parameters n, gnn, and
gn,n1 exhibit a common pattern. If we take the
ratio of a sensitivity in the reference system to the corresponding
sensitivity in the alternative system, we find that the ratio of the
sensitivities is always less than 1. That is,
|
(30)
|
where 
is a positive sum of products of the kinetic orders,
g1n < 0, gj+1,j > 0 for j = 1, 2, ... , n 1, and
|
(31)
|
Thus, the flux in the reference system is less sensitive to
parameter variations, i.e., is more robust than that in the alternative system.
The sensitivities of the steady-state flux with respect to the
parameter g1,n1 exhibit a similar pattern. The
ratio of the sensitivities in this case is given by
|
(32)
|
Although the function of the kinetic orders is different from that
in Eq. 30, the flux in the reference system is again less sensitive to
parameter variations, i.e., is more robust than that in the alternative system.
In contrast, the ratio of the sensitivities with respect to the
parameters n+1 and gn+1,n
exhibits a different pattern,
|
(33)
|
where is a negative sum of products of the kinetic orders.
These parameter sensitivities are related to the last enzyme and
reflect the design for responsiveness to changes in demand for end product.
As the position of the last intermediate that provides feedback
inhibition to the first reaction approaches the beginning of the
pathway, the number of sensitivities that differ between reference and
alternative systems increases. This is so because the number of primed
parameters decreases and a smaller number of conditions for external
equivalence are needed to eliminated the extra degrees of freedom. In
general, if the last intermediate that provides an inhibitory feedback
to the first reaction is Xk for k < n 1, then the sensitivities of the flux to the rate constants k to n+1 and those to the
kinetic orders gij (k i n and i j n) will differ between the
reference and the alternative systems. In most cases, the sensitivities
will be less in the reference system. There are exceptions to this, depending on the length of the pathway and on the last intermediate that provides feedback inhibition to the first step, and, in the case
of n+1 and gn+1,n, the
sensitivities of the reference system will always be greater, for the
reasons we have already mentioned.
Robustness of concentrations
The steady-state concentrations of reference and alternative
systems have different sensitivities with respect to many parameters that define the systems. In some cases, the ratio of the corresponding sensitivities is always <1 or always >1, but, in others, the ratio is
<1 for some values of the parameters and >1 for other values. In the
latter cases, an examination of actual numerical values for the
parameters is critical.
The ratio of sensitivities for the concentration of each intermediate
in the pathway with respect to changes in the kinetic order
g1,n1 is identical to that given in Eq. 32.
Similarly, the ratio of sensitivities for Xn
with respect to changes in the rate constants n or
n+1 is always of the form
|
(34)
|
where p is a different positive sum of products of
kinetic orders for each p, p = n, n + 1, and
|
(35)
|
Thus, the reference system is always less sensitive to changes in
these parameters.
In contrast, the ratio of sensitivities for Xn
with respect to changes in the kinetic orders
gn+1,n, gn,n1, or gnn is always of the form
|
(36)
|
where pq is a different positive sum of products of
the kinetic orders for each gpq. In this case,
the ratio can be >1 or <1. This means that the sensitivity of the
reference system will be greater than the sensitivity of the
alternative system for some values of the parameters and less for
others. Similarly, the ratio of sensitivities for each intermediate
Xi with respect to changes in each parameter can
be 1, depending on values of the parameters.
Again, as the position of the last intermediate that provides feedback
inhibition to the first reaction approaches the beginning of the
pathway, the number of sensitivities that differ between reference and
alternative systems increases. In general, if the last intermediate
that provides an inhibitory feedback to the first reaction is
Xk, then the ratio of sensitivities for each metabolite with respect to changes in the kinetic order
g1k is given by
|
(37)
|
In this equation, 1k is a positive subdeterminant
of the [A] matrix. The ratio of sensitivities for the end product
with respect to changes in each of the parameters common to the two systems also is always 1. Similarly, the ratio of sensitivities for
the last intermediate that feeds back to the first reaction, Xk, with respect to the parameters
k or gkj (k j n) is always <1. Thus, the reference system is always more
robust than the alternative system in these cases. As for the remaining
cases, the sensitivities of the reference system will be greater than the sensitivities of the alternative system for some values of the
parameters and less for others.
Stability
The characteristic equation for Eqs. 1-3 operating near the
steady state can be written as
|
(38)
|
where Fi = Vi0/Xi0 and
aij = gij gi+1,j. Eq. 38 can be expanded into polynomial form
and the Routh conditions for local stability determined. The last two
Routh conditions are critical for stability (Frazer and Duncan,
1929). The last condition is equivalent to the condition
(1)ndet(A) > 0, which is always true for the
systems we are considering (Savageau, 1976, Appendix B).
The two critical Routh conditions for a two-step pathway are
|
(39)
|
and
![R<SUB>2</SUB>=F<SUB>1</SUB>F<SUB>2</SUB>[g<SUB>11</SUB>(g<SUB>22</SUB>−g<SUB>32</SUB>)+g<SUB>21</SUB>(g<SUB>32</SUB>−g<SUB>12</SUB>)]>0.](/content/vol79/issue5/fulltext/2290/img040.gif)
|
(40)
|
Both these conditions are always satisfied for both system A
(g12 < 0) and system B
(g'12 = 0 and
g'11 = g11 + g12g21/(g32 g22) < g11 < 0), so these
systems are always stable. The ratio of the last Routh condition for
the two systems is equal to unity, whereas that for the penultimate
condition is given by
|
(41)
|
Thus, the stability margin is larger for the alternative
system B.
The two critical Routh conditions for a three-step pathway are already
considerably more complex. Whereas the last condition is always
positive, the most critical condition is the penultimate one that can
be positive or negative, depending upon the particular values for the
parameters. The ratio of the last condition for the two systems is
equal to 1; the ratio of the penultimate condition can be >1 or <1,
depending on the values for the parameters. These same conclusions are
obtained for pathways of length four or greater: the ratios cannot be
determined analytically to be >1 or <1, and we must resort to
numerical methods.
Transient time
There is no analytical way to accurately calculate the transient
times of the pathway. This must be done numerically.
Numerical comparisons
Unlike the symbolic analysis performed in the previous section,
using actual numbers for the values of the parameters limits the
absolute generality of the results. However, it does allow us to obtain
general conclusions in a statistical sense. The results described below
have been obtained for pathways of up to seven intermediates. The
trends in these results remain constant throughout all the tested
lengths (i.e., pathways from 2 to 7 intermediates), which suggests that
they will remain so for longer pathways. The use of these numerical
methods allows us not only to study the effects of overall feedback
inhibition, but also to study correlations that exist between systemic
properties and the different parameters of the system.
Response to availability of substrate and demand for end
product
The logarithmic gains in concentrations of the two systems in
response to changes in the initial substrate X0
are identical at each step in the pathway because of the constraints
for external equivalence described in the Methods section. The same is
true for the logarithmic gains in flux. Hence, the responsiveness of concentrations and fluxes in the two systems to changes in the availability of initial substrate is identical numerically as well as analytically.
The logarithmic gain in flux for system A in response to changes in the
demand for end product was shown analytically to be greater than that
for system B. The graph of
L(V, Xn+1)A/L(V, Xn+1)B
versus L(V, Xn+1)A (Fig.
2 A), which is the moving median density of ratios plot introduced in Alves and Savageau (2000c), shows how much greater, on average, the response is
for system A. It also shows a negative correlation between the ratio of
responses and the response of the reference system. This means that, as
L(V, Xn+1)A increases, the ratio
L(V, Xn+1)A/L(V, Xn+1)B
tends to decrease.
View larger version (34K):
[in this window]
[in a new window]
|
FIGURE 2
Typical moving median density of ratios plots for
different magnitudes. The values on the X-axis represent the
moving median of the relevant magnitude in the reference system. The
values on the Y-axis represent the moving median of the
ratio of that magnitude in the reference system to the corresponding
magnitude in the alternative system. (A) Logarithmic gain in
flux in response to changes in demand for the end product,
L(V, Xn+1). (B) Logarithmic gain in
end-product concentration in response to changes in demand for the end
product, L(Xn, Xn+1).
(C) Aggregate sensitivity of the pathway flux,
S(V). (D) Aggregate sensitivity of the
concentration of the last intermediate to feed back on the first
reaction, S(Xk). (E) Aggregate
sensitivity of the concentration of any intermediate in the pathway
before Xk, S(Xi).
(F) Aggregate sensitivity of the concentration of any
intermediate in the pathway after Xk,
S(Xj). (G) Aggregate sensitivity of the
concentration of the end-product, S(Xn).
(H) The penultimate (i.e., n 1st) Routh
criterion; this represents the margin of stability. (I)
Transient time,  in normalized units, is the time the pathway takes
to return within 1% of its steady state following a 15% perturbation
in the steady-state values. Each of these plots is for a specific
pathway length; only the parameter values are changed randomly.
However, because the trends observed for different pathway lengths are
the same, we have only shown a representative case.
|
|
The logarithmic gain in end-product concentration for system A in
response to changes in the demand for end product also was shown
analytically to be smaller than that for system B. The graph of
L(Xn, Xn+1)A/L(Xn, Xn+1)B
versus
L(Xn, Xn+1)A (Fig. 2 B) shows how much smaller, on average, the response
is for system A. It also shows a positive correlation between the ratio
of responses and the response of the reference system.
Robustness
Figure 2 shows typical moving median density of ratios plots for
the aggregate parameter sensitivities of flux and concentrations. The
aggregate parameter sensitivity of the flux V is smaller, on
average, for system A (Fig. 2 C). Assume that
Xk is the last intermediate to feed back on the
first reaction of the pathway. The aggregate parameter sensitivity of
Xk is smaller, on average, for system B (Fig.
2 D). The average difference in aggregate sensitivities for
this metabolite is never larger than a few percent. With regard to the
remaining intermediates, the graphs for Xi (Fig.
2 E) and Xj (Fig. 2 F)
represent typical plots of aggregate parameter sensitivities. In these
cases, we find that random reference systems are less sensitive than
the equivalent alternative systems. The average differences can range
from a few percent to fifty or more percent. The individual parameter
sensitivities of Xn were analytically determined
to be smaller in system A. In the example presented here, the
difference is, on average, just a few percent (Fig. 2 G);
however, depending on the length of the pathway, this
difference can increase to more significant values.
The flux (Fig. 2 C) and concentrations
Xi, i < n, (Fig. 2,
D, E, and F) show a positive
correlation between the ratio of their aggregate sensitivities in the
two systems and the aggregate sensitivity in the reference system when
its value is low. For systems with low sensitivities, system A is, on
average, much less sensitive than system B. For higher values of the
aggregate sensitivities in the reference system, there is no
correlation. In the case of Xk, the ratio is
fairly independent of the values of the aggregate sensitivity in the
reference system.
Stability
The last critical Routh criterion is always the same in
the reference and alternative systems, as has been shown analytically. For a two-step pathway, the margin of stability determined by the
penultimate criterion is always larger in system B. For longer pathways, the margin of stability can be larger in either the reference
or the alternative system, depending on the numerical values of the
parameters. The differences between the two systems with respect to
this penultimate criterion are small (on average less than 2%,
Fig. 2 H), which implies that systems with and
without overall feedback inhibition will have comparable stability margins.
Transient time
Fig. 2 I shows a typical moving median
density of ratio plot for transient time. This plot shows that the
reference system usually responds to perturbations in the steady state
more quickly than the alternative system. For reference systems with a
fast response to changes, the transient times can be, on average, half that of the corresponding alternative systems. For reference systems that are sluggish, the difference is, on average, smaller, though it
still exists.
Effects of parameter values on systemic properties
Rate-constant effects on aggregate sensitivities
Assume that Xk is the last
intermediate to feed back on the first reaction. Plotting the aggregate
sensitivities as a function of j, n j, shows that there is a correlation between each rate constant
j and each of the aggregate sensitivities (Fig.
3 A). For small
j, the correlation is either nonexistent or slightly negative, whereas, for large values, this correlation is positive. As
for the other rate constants, with j < n, there are no
obvious correlations that are general for all the pathway lengths
studied, although, for some lengths, specific correlations are
observed.
View larger version (31K):
[in this window]
[in a new window]
|
FIGURE 3
Typical moving median correlation plots between
different systemic properties and different kinetic parameters of the
reference system. The values on the X-axis represent the
moving median of the relevant kinetic parameter. The values on the
Y-axis represent the moving median of the relevant systemic
property. (A) Aggregate sensitivity of the concentration of
any pathway intermediate Xi versus the
rate-constant parameters n or n+1.
(B) Aggregate sensitivity of the concentration of the end
product Xn versus the kinetic-order parameters
g1n or gi+1,n.
(C) Aggregate sensitivity of the concentration of any
pathway intermediate Xi versus the kinetic-order
parameters g1+1,i or
gn,n1. (D) Aggregate sensitivity of
the concentration of any pathway intermediate Xi
versus the kinetic-order parameter gii.
(E) Aggregate sensitivity of the concentration of the end
product Xn versus the kinetic-order parameter
gn+1,n. (F) Aggregate sensitivity of
the pathway flux V versus the kinetic-order parameter
gn+1,n. (G) Aggregate sensitivity of
the pathway flux V versus the kinetic-order parameter
gn,n1. (H) Transient time versus the kinetic-order parameter g1i.
(I) Transient time versus the kinetic-order parameter
gi+1,i. Each of these plots is for a specific
pathway length; only the parameter values are changed randomly.
However, because the trends observed for different pathway lengths are
the same, we have only shown a representative case.
|
|
Kinetic-order effects on aggregate sensitivities
For Xn, the aggregate sensitivity is
correlated with several parameters. There is a positive correlation
between this sensitivity and g1n. Because
g1n is always negative, this means that the
aggregate sensitivity of Xn,
S(Xn), is usually smaller for high values of overall
feedback inhibition. The same is true for the correlation between
S(Xn) and gin when
i < n (Fig. 3 B). If i = n, there is a negative correlation between this aggregate
sensitivity and g1n. The correlation of the
aggregate sensitivities of the other intermediates with
g1n is usually small or nonexistent. There is a
negative correlation between the aggregate sensitivity of Xi and gi+1,i or
gn,n1 (Fig. 3 C) and a positive
correlation between that of Xi and
gii (Fig. 3 D). Also, the aggregate
sensitivity of each X is negatively correlated with
gn+1,n (Fig. 3 E). These are the
correlations that are generally observed for the aggregate
sensitivities of concentrations, although other individual correlations
can be found for specific intermediates and specific pathway lengths.
The correlations between aggregate sensitivities of flux and the
various kinetic-order parameters are less clear. The correlation with
gn+1,n is positive for low values of
gn+1,n, but it disappears as the value of
gn+1,n increases (Fig. 3 F). The
only other general correlation observed is that between the aggregate
sensitivity of the flux and the kinetic order
gn,n1. This is a negative correlation that
also vanishes as the value of gn,n1 increases.
This can be seen in Fig. 3 G.
Rate-constant and kinetic-order effects on margin of stability
The correlations between a given Routh criterion and the various
parameters depends on which criterion is considered. The results are
pathway length-specific, and no general trend can be found.
Rate-constant and kinetic-order effects on transient time
There is no clear correlation between transient time and the
various rate constants. There are, however, positive correlations between transient time and the kinetic orders
g1i, for i 
1 (Fig. 3 H). There also are negative correlations between
transient time and the kinetic orders gi+1,i for
i > 1 (Fig. 3 I). These were the only
observed correlations with transient time.
Effects of enzyme levels on systemic variables
We have determined the logarithmic gains in flux and
concentrations in response to changes in the level of individual
enzymes. When comparing logarithmic gains in flux and concentrations in the reference and alternative systems, the equivalence conditions will
make all corresponding coefficients identical except the last two. We
also have examined the correlations among the logarithmic gains.
The last two logarithmic
and
TOP
ABSTRACT
INTRODUCTION
