Emergence of self-sustained oscillations in a hypothetical metabolic pathway regulated by a positive feedback loop

Introduction

1.1 Oscillations in chemical and biological systems

In the early 1950s, Russian biochemist Boris Pavlovich Belousov discovered that the oxidation of citric acid, with bromate (BrO3-) and cerium ions as a catalyst, did not progress as expected. Instead of a gradual change from a colorless solution to a pale yellow one, due to the oxidation of Ce3+ to Ce4+, periodic fluctuations in the color of the reaction mixture were observed, indicating that successive interconversions between reduced and oxidized states of reactants were taking place. Belousov was unsuccessful in publishing his results because the reviewers believed that this phenomenon disagreed with the Second Law of Thermodynamics (Kiprijanov, 2016). However, in 1961 Zhabotinsky verified Belousov’s discovery, publishing his results in 1964 (Tyson, 1994). The Belousov-Zhabotinsky reaction is considered a prototype of an oscillatory chemical reaction, being its mechanism well-known (Noyes et al., 1972). Nevertheless, oscillations in chemical systems were discovered long before. As early as the 17th century, Robert Boyle observed the intermittent ignition of phosphorus (Agreda & Barragán, 1998), and references from the 19th century mention the mercury beating heart, an electrochemical-based mechanical oscillator (Martín & Martín, 2010). In 1910 Lotka proposed a hypothetical reaction model that displays damped oscillations (Lotka, 1910).

Later on, in 1920, the same author showed that chemical systems could adopt sustained oscillatory behavior under certain circumstances (Lotka, 1920). However, it was not until the late 1960s when other researchers became interested in these types of reactions (Goldbeter, 2018). In 1956 Prigogine and Balescu demonstrated that sustained oscillations could occur in open chemical systems with nonlinear kinetics and far from equilibrium (Prigogine & Balescu, 1956; Goldbeter, 1996). In the 1960s, Goodwin studied negative feedback as a source of oscillations and the quasiperiodic and subharmonic behavior resulting from the interaction between two oscillators. Goodwin ventured the importance of these phenomena for cellular physiology (Gonze & Ruoff, 2020). In 1967 Prigogine and Nicolis demonstrated that sustained oscillations do not contradict any thermodynamic principle and showed that the instability of these chemical oscillators is due to positive (PFBLs) or negative feedback loops (NFBLs) (Prigogine & Nicolis, 1967; Agreda & Barragán, 1998).

In 1968 Lefever, Nicolis and Prigogine presented an autocatalytic theoretical model known as Brusselator (an acronym for Brussels and oscillator) that has been very useful for the understanding of oscillatory behaviors in chemical systems (Lefever et al., 1967; Goldbeter, 1996; Agreda & Barragán, 1998). In contrast to the scarcity of oscillatory behavior in chemical systems, rhythmicity is inherent to life (Li & Yang, 2018). Although periodic phenomena are not exclusive to biological systems, this behavior is much more frequent in living beings than in purely chemical systems (Goldbeter, 1996; Goldbeter, 2018). The circadian rhythm of sleep and wakefulness is a classic example of rhythmicity in which both frequency and amplitude of a circadian oscillator must keep robust 24-hour rhythms in the presence of external fluctuations (Kang & Cho, 2017). Although oscillations are necessary for many biological functions that require precise control of time, biological systems are inherently noisy, affecting its oscillations’ accuracy.

However, these systems achieve robustness by dissipating a critical amount of free energy, which results in greater irreversibility and accuracy of the oscillations (Cao et al., 2015; Fei et al., 2018). The relationship between the amount of energy dissipated and the precision of the oscillations period is a particularly topical question that has been the subject of recent research (del Junco and Vaikuntanathan, 2020). Heart and respiratory rates, fluctuations in hormone levels, electrically excitable cells like neurons, or muscle cells are additional examples of rhythmicity. Physiological oscillations can be a consequence of the central nervous system’s activity, although in certain instances, they are generated in situ. Thus, the heart rate is due to the sinoatrial node’s pacemaker cells, located in the right atrium, which have an unstable resting membrane potential.

Rhythmically and spontaneously, these cells generate an action potential that is transmitted with the necessary delay to coordinate the atria and ventricles’ function throughout the myocardium whose cells are electrically coupled. The frequency of the oscillations can become high; for example, in those insects equipped with asynchronous flying muscles, as Diptera, the wings can reach, in certain species, a frequency of 1000 Hz (Ruppert & Barnes, 1996). Oscillatory behavior is also observed in non-excitable cells, such as intracellular Ca2+ fluctuations, enzymatic reactions, or the secretion of many hormones released from the endocrine glands in a pulsatile manner rather than continuously (Goldbeter, 2018). Additional examples of oscillatory behaviors in non-excitable cells are the nuclear factor kappa B (NF-κB), and the proteins p53 and p38 (Li & Yang, 2018). The latter is a mitogen-activated protein kinase (MAPK), which controls an inflammatory response; oscillations in its activity have been observed for more than 8 hours under stimulation by interleukin-1β (Tomida et al., 2015).

The tumor suppressor protein p53, regulated through a negative feedback loop (NFBL) with the mdm2 protein, also exhibits an oscillatory dynamic. This mechanism allows the cell to repair its DNA through the action of p53 without excessive exposure to said protein, which could induce the apoptotic cell-death (Bar-Or et al., 2000). The oscillations of p53 are observed in cell lines of different species, although with different frequencies. For example, in rodent cells, these oscillations are faster than those observed in cell lines derived from humans, monkeys, and dogs (Stewart-Ornstein et al., 2018). Nuclear factor-kappa B (NF-κB), a protein complex that plays a crucial role in inflammation and cancer, exhibits damped oscillations synchronized with external disturbances (Zambrano et al., 2016).

1.2 Biochemical oscillations

The study of oscillating biochemical systems began in 1955 after observing changes in the concentrations of 3-phosphoglycerate and ribulose-1,5-diphosphate during the dark phase of photosynthesis (Chance et al., 1973). Oscillations in glycolysis were observed in 1957, while continuous glycolytic oscillations in cell-free yeast systems were demonstrated in 1966. (Chance et al., 1973; Agreda & Barragán, 1998; Robles & Barragán, 2006). Oscillations of cAMP (Li & Goldbeter, 1990; Goldbeter, 2018), intracellular Ca2+ (Dupont & Goldbeter, 1992), and glycolysis metabolites (Goldbeter, 1996) are among the best-known examples of periodic behavior at the biochemical level. The role of oscillations can be diverse and have high functional importance; abnormalities in them can cause disease, including cancer (Li & Yang, 2018). It has been suggested that they may be at the origin of some circadian rhythms and other processes such as cell division (Rodríguez et al., 1998). Thus, in Dyctiostelium discodeium amoebas, concentration’s periodic oscillations are related to the aggregation of cells that respond chemotactically to cAMP pulses (Goldbeter 1996; Goldbeter, 2018).

Periodic and chaotic behaviors are observed in metabolic pathways driven by matter flows and controlled by positive or negative feedback loops (PFBLs or NFBLs) (Goldbeter, 1996). In the context of biochemical pathways, feedback loops (FBLs) consist of at least two metabolites that regulate their activity reciprocally by activation or inhibition. The complex behavior of biological systems results from the coexistence of multiple of these FBLs (Dalchau et al., 2018).

Glycolysis is considered the central route of glucose catabolism and provides cells a significant portion of energy and pyruvate that can continue the subsequent oxidative degradation. It is a 10-steps metabolic pathway that transforms glucose and other hexoses into two pyruvate molecules, obtaining two molecules of ATP and NADH. This metabolic pathway begins and ends with irreversible steps: the conversion of glucose to glucose-6-phosphate by hexokinase, consuming ATP (step 1), and the conversion of phosphoenolpyruvate to pyruvate by pyruvate kinase, producing ATP (step 10). There is a third irreversible step: the transformation of fructose-6-phosphate to fructose-1,6-bisphosphate by phosphofructokinase (PFK). These three irreversible steps suppose a greater dissipation of energy, driving glycolysis away from thermodynamic equilibrium (Robles & Barragán, 2006).

Metabolic oscillations in this pathway have been observed in tumor cells, muscle extracts, fibroblasts, and yeast, the latter being the most studied system (Sel’kov, 1968; Robles & Barragán, 2006). These oscillations emerge in glycolysis when glucose-6-phosphate or fructose-6-phosphate is supplied, but not for the addition of fructose-1,6-diphosphate. This observation suggests that the enzyme phosphofructokinase (PFK) is the oscillophore responsible for the periodic behavior of the pathway, which can also be demonstrated by adding inhibitors or activators of the enzyme. Specifically, the addition of ADP (an activator of PFK) affected the frequency and amplitude of the oscillations of glycolytic metabolites; this is an example of PFBL.

Although PFBLs exerted by a product of a metabolic pathway on an enzyme are not as frequent as product inhibition, autocatalysis has an essential role in biochemical oscillations (Robles & Barragán, 2006; Tyson, 2002). Nonetheless, the importance of PFK in the oscillations observed in this metabolic pathway has also been relativized, proposing that the characteristics of the oscillations also depend on other enzymes of the pathway (Reijenga et al., 2005). Alternatives to the oscillophore have been proposed to explain this glycolytic behavior through the autocatalytic character of the stoichiometry of the pathway, which acts as a PFBL responsible for the oscillations (Gustavsson et al., 2014).

The theoretical approach to these periodic biochemical phenomena through models is useful in understanding their physiological roles and underlying mechanisms. An example of the latter is the Oregonator, an autocatalytic model that reproduces the oscillatory behavior of the Belousov-Zabotinsky reaction (Field & Noyes, 1974).

This article presents a hypothetical metabolic pathway with a PFBL. The system is kept away from thermodynamic equilibrium by constant input and output fluxes. The positive feedback consists of the activation of the first enzyme of the pathway by its product. Our main target is to provide a simple kinetic model in which, intrinsically and under certain conditions, sustained oscillatory behavior arises. Such a model could be a starting point for building modified or more complex and realistic models for educational and research purposes.

2 Model and kinetic equations

Figure 1 shows the sequence of enzymatic reactions of a hypothetical linear metabolic pathway consisting of four steps and six metabolites: S0, S1, S2, S3, S4, and S5. The concentration of S0 is constant and it is also converted at a constant rate by a non-enzymatic irreversible step in S1, representing a non-rhythmic input of the substrate into the system. The rest of the metabolites are sequentially transformed from S1 to S5 through the subsequent action of the enzymes E₁, E₂, E₃, and E₄.

$$\begin{gathered} {\mathrm{S}}_{\mathrm{o}}\stackrel{{\mathrm{k}}_{\mathrm{in}}}{\to }{\mathrm{S}}_{1 }+ {\mathrm{E}}_1\underset{{\mathrm{k}}_{-1}}{\overset{{\mathrm{k}}_1}{\rightleftharpoons }}{\mathrm{E}}_1{\mathrm{S}}_1\stackrel{{\mathrm{k}}_2}{\to }{\mathrm{E}}_1+{\mathrm{S}}_2 \\ {\mathrm{S}}_{2 }+ {\mathrm{E}}_2\underset{{\mathrm{k}}_{-3}}{\overset{{\mathrm{k}}_3}{\rightleftharpoons }}{\mathrm{E}}_2{\mathrm{S}}_2\stackrel{{\mathrm{k}}_2}{\to }{\mathrm{E}}_2+{\mathrm{S}}_3 \\ {\mathrm{S}}_{3 }+ {\mathrm{E}}_3\underset{{\mathrm{k}}_{-5}}{\overset{{\mathrm{k}}_5}{\rightleftharpoons }}{\mathrm{E}}_3{\mathrm{S}}_3\stackrel{{\mathrm{k}}_6}{\to }{\mathrm{E}}_3+{\mathrm{S}}_4 \\ {\mathrm{S}}_{4 }+ {\mathrm{E}}_4\underset{{\mathrm{k}}_{-7}}{\overset{{\mathrm{k}}_7}{\rightleftharpoons }}{\mathrm{E}}_4{\mathrm{S}}_4\stackrel{{\mathrm{k}}_8}{\to }{\mathrm{E}}_4+{\mathrm{S}}_5\stackrel{{\mathrm{k}}_{\mathrm{out}}}{\to } \end{gathered}$$

Figure 1. Sequential reactions that make up a linear metabolic pathway

As described in Figure 2, each enzyme reversible binds to its substrate in these four enzymatic reactions, forming a substrate-enzyme complex (E_x S_x ) that converts into the corresponding product. The kinetic constants for the reversible formation of the substrate enzyme complex are given by k_n and k_-n. For each step of the pathway, the substrate enzyme complex (Ex Sx ) yields its corresponding product (S_x+1) with a kinetic constant denoted by k_{n + 1} and k_(n+1)-1 for the reverse reaction. When k_(n+1)-1= 0, the conversion of E_x S_x into E_x and S_x+1 is irreversible. 

${\mathrm{E}}_{\mathrm{x}}+{\mathrm{S}}_{\mathrm{x} }\underset{{\mathrm{k}}_{-\mathrm{n}}}{\overset{{\mathrm{k}}_{\mathrm{n}}}{\rightleftharpoons }}{\mathrm{E}}_{\mathrm{x}}{\mathrm{S}}_{\mathrm{x}}\underset{{\mathrm{k}}_{(\mathrm{n}+1)-1}}{\overset{{\mathrm{k}}_{\mathrm{n}+1}}{\rightleftharpoons }}{\mathrm{E}}_{\mathrm{x} }+ {\mathrm{S}}_{\mathrm{x}+1}$

Figure 2. Enzyme reversibly binds its substrate forming a substrate-enzyme complex

Note: Each enzyme of the pathway reversibly binds its substrate to form an enzyme-substrate complex (E^xS^x) with kⁿ and k^-n as kinetic constants of the direct and reverse reaction. The E^xS^x complexes evolve to form their corresponding products (S^x+1) with a constant rate of ^kn+1.

The dynamic of the enzymatic conversions showed in Figure 1 can be described by the following differential equations:

$1) \frac{\mathrm{d}\left[{\mathrm{S}}_1\right]}{\mathrm{dt}}= {\mathrm{k}}_{\mathrm{in}}\left[\mathrm{S}0\right]-{\mathrm{k}}_1\left(1+{\left[{\mathrm{S}}_2\right]}^3\right)\left[{\mathrm{E}}_1\right]\left[{\mathrm{S}}_1\right]+{\mathrm{k}}_{-1}\left[{\mathrm{E}}_1{\mathrm{S}}_1\right]$

$2) \frac{\mathrm{d}\left[{\mathrm{S}}_2\right]}{\mathrm{dt}}= {\mathrm{k}}_2\left(1+{\left[{\mathrm{S}}_2\right]}^3\right)\left[{\mathrm{E}}_1{\mathrm{S}}_1\right]-{\mathrm{k}}_3\left[{\mathrm{E}}_2\right]\left[{\mathrm{S}}_2\right]+{\mathrm{k}}_{-3}\left[{\mathrm{E}}_2{\mathrm{S}}_2\right]$

$3) \frac{\mathrm{d}\left[{\mathrm{S}}_3\right]}{\mathrm{dt}}= {\mathrm{k}}_4\left[{\mathrm{E}}_2{\mathrm{S}}_2\right]-{\mathrm{k}}_5\left[{\mathrm{E}}_3{\mathrm{S}}_3\right]+{\mathrm{k}}_{-5}\left[{\mathrm{E}}_3{\mathrm{S}}_3\right]$

$4) \frac{\mathrm{d}\left[{\mathrm{S}}_4\right]}{\mathrm{dt}}= {\mathrm{k}}_6\left[{\mathrm{E}}_3{\mathrm{S}}_3\right]-{\mathrm{k}}_7\left[{\mathrm{E}}_4{\mathrm{S}}_4\right]+{\mathrm{k}}_{-7}\left[{\mathrm{E}}_4{\mathrm{S}}_4\right]$

$5) \frac{\mathrm{d}\left[{\mathrm{S}}_5\right]}{\mathrm{dt}}= {\mathrm{k}}_8\left[{\mathrm{E}}_4{\mathrm{S}}_4\right]-{\mathrm{k}}_{\mathrm{out}}\left[{\mathrm{S}}_5\right]$

As shown in Figure 1, the E1-catalyzed transformation of S₁ to S₂ simulates an activation per product step. The formation reaction rate for the E₁S₁ complex and its conversion into E₁ and S₂ have become strongly dependent on the concentration of S₂ by the inclusion of the algebraic expression 1+ [S₂]³ in equations 1 and 2. This represents an activation of the E₁ enzyme by the substrate S₂, increasing both the affinity of the enzyme E1 to the substrate S1 and the conversion rate of this complex to S2 and E1. The inclusion of the binomial 1+ [S₂]³ in the differential equations is responsible for the system’s non-linearity. According to the Arrhenius equation, the kinetic constants depend not only on the temperature but also on the presence of catalyst; therefore, it can be considered that the expression (1+S₂³) modifies the kinetic constants as follows: k₁’ = k₁ (1+S₂³) and k₂’ = k₂ (1+S₂³).

Since in living cells, the inlet flow is considered a key parameter to trigger the oscillations (Danø et al., 1999); the model›s behavior was evaluated for concentration values of S₀ between 0 and 12 arbitrary units (a.u.). In all experiments, the inlet flux was such that the concentration of S₂ (activator of E₁) was lower than the concentration of E₁ ([S₂] < [E₁]). Table 1 shows the values of the kinetic parameters, concentration of enzymes, and initial substrates concentration for each of the pathway steps.

Python version 3.7.4 with Spyder and Anaconda as IDE was used for numerical integration of differential equations. The simulations were performed, implementing the kinetic parameters and concentration values for metabolites and enzymes shown in Table 1. Python script used for the experiments is provided in the supplementary materials.

Results

Using the values specified in Table 1 and for [So] equal to 8 a.u., the system reaches the steady state (Figure 3) with saturation percentages for enzymes not exceeding 40%. The steady state concentrations of metabolites and enzyme saturation percentages are shown in Table 2.

However, when [S₀] = 10 a.u., the system does not reach a steady state but is installed in a sustained oscillation, as displayed in Figure 4. Periodic oscillations with T=1.775 a.u. are observed in all metabolites of the pathway. The transformation of S₁ to S₂ speeds up as S₂ increases until S₁ depletion (Figure 1), limiting S₂ production to a value that allows recovery of S₁ levels, this being a classic oscillation generator mechanism (Tyson, 2002).

The oscillatory ranges of the concentration of metabolites are shown in Table 3.

Figure 5 illustrates the phase plane plot for S₂ against S₃, observing a limit cycle that is the proper attractor for the oscillatory behavior.

The transition from the steady state to the sustained oscillation takes place through a Hopf bifurcation (Gustavsson et al., 2014; Beta & Kruse, 2017). For a critical value of S₀ (approximately 9.8 a.u.), the steady state becomes unstable, and the system evolves to an oscillatory behavior. When S₀=9.7 a.u., near the bifurcation point, the system shows damped oscillations, which eventually reaches the steady state as shown in Figure 6. Figure 7 shows the phase plane for this experiment; a stable focus is attained. Figure 8 displays a close-up view of the said attractor.

The system’s evolution from a fixed-point attractor towards the limit cycle has also been studied (Figure 9). As previously explained, when the concentration of S₀ is equal to 9.7 a.u. the limit cycle is not achieved, resulting in a stable focus (Figure 9A). However, starting from a concentration of S₀ equal to 9.8 a.u., a stable limit cycle starts taking shape as a result of the self-sustained oscillations, as described by Saha & Gangopadhyay (Saha & Gangopadhyay, 2017). This tendency is confirmed for concentration values of S₀ equal to 9.9 and 10 a.u., for which their corresponding phase planes are presented in Figure 9, labeled as (C) and (D). The system returns to a steady state when the concentration of S₀ is equal or higher than 11.75 (a.u.), as demonstrated in Figure 10. It should be clarified that the emergence of oscillatory behavior also depends on the concentrations of the enzymes. For example, changing the concentrations of E₁, E₂, E₃, and E₄ to 15 a.u. and [S₀] equal to 10 a.u., the system does not oscillate. Instead, it reaches a steady state, thus making the oscillatory window even narrower. Contrary to what happens within the oscillatory interval, the concentration of S₂ did not drop enough to allow a significant increase in S₁, which remains at low concentration, being this the probable reason for the lack of oscillations.