Royal Society Publishing

Collapse of single stable states via a fractal attraction basin: analysis of a representative metabolic network

Junli Liu, John W Crawford, Konstantinos I Leontiou


The impact of external forcing on an enzymatic reaction system with a single finite stable state is investigated. External forcing impacts on the system in two distinct ways: firstly, the reaction system undergoes a series of discontinuous changes in dynamical state. Secondly, a critical level of forcing exists, beyond which all finite states become unstable. It is shown that the results stem from the conditions for global stability of the system. Competition between the attractor for stable states and the unbounded states leads to a loss of integrity and the fractal fragmentation of the attraction basin for the finite state. The consequences of a fractal basin in this context are profound. Initial states which are infinitesimally close diverge to a finite and an unbounded state where only the finite state is consistent with biological functionality. Furthermore, above a critical forcing amplitude, the system does not converge to a finite state from any initial state, implying that there is no configuration of metabolite concentrations that is consistent with sustained evolution of the system. These results point to opportunities for constraining uncertainty in cell networks where nonlinear saturating kinetics form an important component.


1. Introduction

The study of the rich dynamical behaviour that results from nonlinear reaction kinetics has yielded new insights not only into the interplay of reactions but also in dynamical systems theory itself (Scott 1991; Goldbeter 1996). Complex dynamical behaviour, such as periodic and chaotic oscillations, as well as the coexistence of multiple stable attractors has been observed in both chemical and biochemical reactions (Scott 1991; Goldbeter 1996). Where multiple stable attractors coexist, each attractor has its own basin of attraction and initial conditions determine which stable state is selected. In some instances, the boundary separating the basin of attraction may become fractal and this leads to a particularly sensitive dependence of state selection on initial conditions (Decroly & Goldbeter 1984; Liu & Scott 1991; Liu et al. 1995; Wang et al. 2003). Ultimately, which of the coexisting stable states is selected is determined by the competition of those stable attractors. Once a state is selected, the system will reach that state asymptotically. This work investigates the loss of integrity of an attraction basin in an enzymatic reaction system with a stable attractor (Liu 1999a) that does not coexist with any other stable attractors. Therefore, biologically, there is only one stable state that can be investigated in the laboratory. The behaviour of this system is analysed in terms of the competition between the stable attractor and unbounded solutions.

In the context of biological pathways, an enzyme-catalysed reaction is embedded in a pathway, taking product molecules of the preceding reaction step and supplying substrate for the subsequent step (Stryer 1997). Enzyme molecules bind with substrate molecules to form an intermediate complex and this step reduces the activation energy of the reaction. Subsequently, the intermediate complex releases product molecules and the enzyme becomes available for catalysis again. Finally, the freed enzyme molecules bind substrate molecules again and the process repeats. Enzyme catalysis can be inhibited or activated by compounds which are in themselves reaction products and the consequent network of feedback and feed-forward reactions is the basis of biological functioning in cells. In a pathway, such a process settles onto a stable, though possibly time-dependant, state. Although metabolic networks can be described in terms of flux control (Fell 1997), spatio-temporal behaviour (Hess 1997) and energy utility (Ross & Schell 1987), their behaviour ultimately arises from the formation of stable states (Liu et al. 1997).

Enzymatic kinetics are usually derived from traditional mass-action kinetics together with simplifying assumptions such as the existence of a quasi-steady state (Segel & Selmrod 1989; Stolerius et al. 2004). The essential feature of the resulting kinetics is a nonlinear dependence of reaction rate on metabolite concentration taking the form of saturation kinetics. This feature leads to the interesting possibility that enzymatically catalysed networks can lose ‘coordination’. By definition, coordination is referred to as the average flux into each species pool being equal to its average output flux (Liu 1999a,b; Liu & Crawford 2000) for a steady state or a time-dependant state. Under mass balance, the coordination of enzymatic reactions may or may not be maintained, depending on the specific relationship among kinetic parameters (Liu & Crawford 2000). The runaway accumulation of metabolites is only prevented if system coordination is maintained and a biochemical system can form stable states with finite metabolites concentrations (Liu 1999a,b; Liu & Crawford 2000). Recently, kinetic modelling for actual biochemical systems has shown that the establishment of flux balances for a steady state requires specific constraints on kinetic parameters. In order to obtain a stable finite state based on the parameters in literature, many of the parameters need to be adjusted (Teusink et al. 2000; Rohwer & Botha 2001; Aon & Cortassa 2002). Even where coordination is guaranteed in an unperturbed system, the addition of an external forcing may result in a subsequent loss of coordination (Liu & Crawford 2000).

In general, all enzymatic reactions in a living system are subject to a fluctuating input. These fluctuations may stem from varying environmental conditions or from the interaction between reactions. For example, plant photosynthesis is subject to light intensity fluctuations and the processes of acquiring carbon resources are time-dependant. Moreover, since a particular reaction is always embedded in a pathway, its input is the output of preceding reactions. By taking into account the complex manner in the network of enzymatic reactions in vivo (Stryer 1997), any reaction is forced by other reactions. For example, oscillations generated in glycolysis can force the pentose phosphate pathway and all enzymatic reactions in that pathway may be subject to a periodic input.

2. Model and the dynamics of an unforced system

In biochemistry, models derived from glycolysis have been intensively studied and represent a class of systems for understanding nonlinear dynamics in metabolic systems (Goldbeter 1996). Recently, a model featuring the conversion of a substrate, S, into a product, P, and into a branched sink has been investigated (Liu 1999a) and is shown in figure 1. In this paper, we use this simple system as a representative example to study the mechanism by which a finite stable state collapses under loss of system coordination.

Figure 1

Reaction scheme. The reaction system includes two molecular species, S and P, and three enzymes, E1, E2 and E3. Enzyme E2 is positively regulated by P.

Biologically, figure 1 represents a part of glycolysis in cells and the underlying assumptions of the model previously fully described by Liu (1999a). Kinetics of all enzymes in figure 1 can be derived based on the basic mass-action law and quasi-steady state assumptions. Recent studies have shown that quasi-steady state assumptions can be generally applied for non-isolated enzymatic reactions described by Michaelis–Menten formalism (Stolerius et al. 2004). Throughout the paper, we assume that kinetics of enzymatic reactions using quasi-steady state assumptions are valid. The governing mass-balance equations are described as follows.Embedded Image(2.1)

Embedded Image(2.2)withEmbedded Image(2.3)Embedded Image(2.4)Embedded Image(2.5)

The biological significances of the parameters were fully described previously (Liu 1999a). In the following, we take Vinput as the bifurcation parameter and fix other parameters (Embedded Image; Embedded Image; Embedded Image; kS=1.0; km2=1.0; kP=1.0; km3=1.0; L=1.0×106).

Following method previously introduced by Davidson & Liu (2002), the global bifurcation diagram with respect to bifurcation parameter Vinput is obtained; it is schematically described in figure 2.

Figure 2

Schematic diagram describing the dynamics of unforced system. There are five regions in terms of the difference in dynamical patterns (see text for details).

There are five different regions in figure 2. (A) A unique finite stable steady state that is globally stable; (B) a unique finite stable steady state that coexists with unbounded solutions; (C) a unique small-amplitude limit cycle coexists with unbounded solutions; (D) a unique large-amplitude limit cycle that coexists with unbounded solutions and (E) unbounded solutions only. The emergence of a large-amplitude limit cycle from a small-amplitude limit cycle through a small increase in Vinput is due to a Canard explosion as exhibited in other systems (Dumortier & Roussarie 1996; Freire et al. 1999; Davidson & Liu 2002). Figure 3 shows the transition from a small- to a large-amplitude limit cycle for two neighbouring values of Vinput.

Figure 3

Comparison of oscillations for two neighbouring input rates: (a) Vinput =2.0580, (b) Vinput =2.0581. In region C (Vinput =2.0580), small-amplitude oscillations emerge. In region D (Vinput =2.0581), large-amplitude oscillations emerge. The transition from region C to region D is due to canard explosion (Davidson & Liu 2002). Initial values: [S]0=1.0; [P]0=1.0.

3. Dynamics and basin of attraction under forcing

In general, the evolution of biochemical systems takes place in the context of fluctuating and changing environments and these fluctuations can take different forms. Here, we assume that the external forcing is represented by the following equation:Embedded Image(3.1)where ϵ and T are the forcing amplitude and period, respectively. In order to guarantee Vinput≥0, ϵ is limited to the range of 0 to 1. We note that the average of Vinput over an exact period is equal to V0, implying that the forcing does not change the net input rate of the system. In this work, we fix V0 to 2.0 so that the system settles into a stable steady state if it is not forced (see figure 2).

(a) Dynamics

In a similar manner to other nonlinear reaction systems (Scott 1991; Goldbeter 1996), the system under forcing supports various periodic and aperiodic oscillations, depending on the parameter values and the forcing period and amplitude. However, the system under study supports additional features, namely, that forcing may lead to unbounded solutions and the consequent loss of system coordination. Figure 4 shows two typical examples for a fixed initial value ([S]0=1.0 and [P]0=1.0). When the forcing period, T=250 s (figure 4a), various periodic and aperiodic patterns emerge. When the forcing amplitude, ϵ, is increased from 0 to 1, no unbounded solutions emerge. However, for T=1500 s (figure 4b), various periodic and aperiodic patterns emerge only if ϵ<0.775 32. When ϵ≥0.775 32, all finite states collapse, leading to unbounded solutions. In general, the system can process external signals across a broad range in a period, provided they have small-amplitude–large-amplitude signals may lead to loss of coordination. Similarly, the system can withstand large-amplitude signals provided the period is short. However, large-amplitude signals with long periods will lead to unbounded solutions. Figure 5 summarizes the effects of forcing amplitude and period on the formation of the stable finite states.

Figure 4

Dependence of dynamics on forcing amplitude for two forcing amplitudes: (a) T=250 s and (b) T=1500 s. In (b), only unbounded solutions exist when ϵ≥0.775 32. Initial values: [S]0=1.0; [P]0=1.0.

Figure 5

Dependence of finite stable-state formation on both forcing amplitude and forcing period. Initial values: [S]0=1.0 and [P]0=1.0.

(b) Basin of attraction

It is clear that, for a fixed initial value, forcing may lead to the collapse of finite states. Biologically, different initial values represent the configuration of the system in the absence of forcing in terms of the different concentrations of metabolites. This configuration of metabolite concentrations can be subject to modification through, for example, genetic engineering. Perhaps more fundamentally important, this configuration might not be known for a particular system under study and may have to be inferred subject to constraints. Clearly, an important constraint is that system coordination must be maintained when the system is subject to perturbation. Therefore, an important question is how different configurations of the system in the absence of forcing could affect the stability of the system, subject to perturbation against loss of coordination. Therefore, we investigated how initial values affect the maintenance of finite states. For the simplicity and without loss of generality, in the following, we fix forcing amplitude to 1500 s.

As the forcing amplitude ϵ increases, there are two distinct regimes that characterize the impact on the size of the finite state attraction basin, as shown in figure 6. For ϵ increasing up to 0.775 31, the size of basin decreases approximately linearly. However, as ϵ increases from 0.775 31 to 0.775 32, the size of the basin decreases abruptly, leading to the loss of all stable finite states. The system does not reach any finite stable state from any initial values when ϵ≥0.775 32.

Figure 6

Dependence of attraction basin of finite stable states on forcing amplitude for T=1500 s. The basin is calculated as follows: the initial value of both S and P are set at zero. Then the initial values of either S or P are increased by 1 each time. If the system develops to a stable finite state, the initial value is recorded. If the system develops to unbounded solutions, the initial value is discarded. The number of the records is the area of attraction basin of finite stable states.

Numerical calculations reveal that the nature of the attraction basin changes qualitatively in the amplitude range of 0.775 31≤ϵ≤ 0.775 32. When ϵ increases up to 0.775 31, the basin of attraction of finite states shrinks gradually but all initial values leading to finite stable states are separated from those that do not by a regular boundary. Therefore, the ‘integrity’ of the basin is maintained. However, as ϵ is gradually increased from 0.775 31 to 0.775 32, the basin of attraction loses its integrity, and the boundary between the finite state basin and the infinite state basin develops a structure which can be approximated by fractal geometry. Figure 7 shows an example of basin of attraction for ϵ=0.775 315 400 1.

Figure 7

Example of fractal boundaries of attraction basin for finite stable state. The stable state does not coexist with any other stable states but with unbounded solutions. Black area leads to finite stable oscillatory states while the corresponding area is left blank for initial values leading to unbounded solutions. The self-similar fine structures at all scales of accuracy reveal the nature of fractal basin boundaries.

Figure 7 shows the attraction basin at different scales of resolution and the self-similar structures of the attraction basin at all scales of accuracy reveal the fractal nature of the basin boundary. For a fractal basin, we would expect the density (ρ(r)) of initial values leading to a finite stable state to depend on the scale of measurement according to the power-law formEmbedded Image(3.2)where r is the scale of measurement and D is the fractal (Hausdorff) dimension. Figure 8 shows a plot of density versus the scale of measurement for the basin shown in figure 7, which indicates the expected linear form (on a double logarithmic plot) predicted by equation (3.2). A linear regression gives an estimate for the fractal dimension of the basin equal to 1.95±0.01 by using a standard statistical package (S-plus). The consequences of this fractal nature of basin boundary are dramatic. Two arbitrarily close initial values may lead to a finite stable state that has biological functionality (figure 9b) and an unbounded solution that cannot perform any function so that the system is biologically unviable (figure 9a). In figure 9b, the concentration of product, [P], is always positive, although it can be low (to four decimal places) at specific stages of the time-dependant response. This reflects that under forcing, the concentration level of product may vary significantly with time.

Figure 8

Dependence of the density of initial values leading to a finite stable state on the scale of measurement (see text for details).

Figure 9

For fractal boundaries in figure 7, two neighbouring initial values may lead to completely different developments of the system. (a) [S]0=74.49 and [P]0=55.41. (b) [S]0=74.49 and [P]0=55.42.

It is clear that external forcing has two distinct impacts on the enzymatic system. Firstly, the system discontinuously changes its dynamical state as the forcing amplitude continuously changes, leading to a range of periodic and aperiodic patterns (figure 4). Secondly, all finite stable states collapse when the forcing amplitude reaches a critical value (figure 6) due to the emergence of a fractal boundary in the attraction basin. Numerical simulations reveal that the collapse of all finite states is not directly associated with the emergence of specific dynamical patterns. The forcing amplitude and period are the key factors determining the collapse of all finite states (figure 6).

4. Concluding remarks

Analysis of a representative metabolic network (figure 1) reveals that, although the system only supports one stable finite state, the stable state may collapse via a fractal basin of attraction under external forcing. When the fractal basin emerges, the evolution of the reaction system is exquisitely sensitive to initial values and two neighbouring initial values can result in divergent evolution: one where the system can implement biological functions and the other where the system is not biologically viable. The phenomenon can be interpreted in terms of the global stability of the system—namely, the coexistence of finite stable states with unbounded solutions.

Biologically, the attraction basin of finite stable states is the range of metabolite concentrations that maintain the stable evolution of the system. External forcing describes environmental fluctuations or stresses. In equation (3.1) the forcing amplitude and period represent the intensity and frequency of environmental fluctuations. Figure 6 clearly shows that, when the forcing amplitude increases up to a critical value (0.775 31), the range of metabolite concentrations for maintaining the stable finite state decreases slowly, implying that the system is robust under the forcing in the sense that large fluctuations in metabolite level can be accommodated. However, when the forcing amplitude increases above a critical value, the range of metabolite concentrations dramatically decreases and only a few selected metabolite concentrations lead to the evolution of stable states. As the forcing amplitude is increased further, no stable states are supported for any metabolite concentrations. Therefore, for any particular configuration of enzyme and metabolite concentrations, the system may only be resistant to a certain level of external forcing. The emergence of a fractal basin of attraction for a unique stable state (in the sense that the state does not coexist with other stable states) has important implications for biological functionality and, therefore, for our understanding of cell networks. There are many approaches to modelling cell networks which vary in the level of detail in the description of the interactions. Inevitably, the most appropriate method is a compromise between realism and the availability of data, with the latter being the ultimate limitation in practice. Therefore, methods which can constrain the uncertainty in parameterization or network configuration are invaluable. The prediction of a fractal basin of attraction and the loss of stable finite states for particular configurations of the network provides theoretical evidence for inconsistent parameterization or network configuration and, thus, provides new constraints which inform our understanding of these kinds of systems. Significantly, these constraints arise from a full nonlinear treatment and cannot be derived from simpler approaches which aim to ignore the inherent nonlinearities of enzymatic systems. This work shows that the special features of nonlinear enzymatic networks have both dynamical and biological implications.

The reaction scheme employed in this work is derived from glycolysis and only includes one positive feedback. Previous studies (Liu 1999b, 2000) have showed that more complicated reaction schemes such as the interplay of positive and negative regulation have similar dynamical features to those in figure 1—namely, finite stable states may collapse. It is expected that the conclusions drawn in this work can be applied to other more complicated systems.

The coexistence of stable states has been found in many physical systems (Thompson 1996). It has been shown (Thompson & Soliman 1990; Sommerer & Ott 1993; Sweet et al. 1999) that in many cases the basin of attraction of a state may change its measure and topology when constraints on the system change. In particular, the basin may become fractal, leading to ‘loss of integrity’ of the state. Under these conditions, the final state of the system depends on the initial state. The distinct feature of the enzymatic reaction system studied in this work is that it has only one finite stable state.


J.L. is supported by the Scottish Executive Environment and Rural Affairs Department (SEERAD).


    • Received May 11, 2004.
    • Accepted October 29, 2004.


View Abstract