Global Hopf bifurcation in networks with fast feedback cycles

TL;DR

This paper identifies structural conditions in networks with feedback cycles that lead to global Hopf bifurcation, resulting in sustained oscillations across various biological and chemical systems.

Contribution

It establishes how monotone feedback cycles in the linearization cause global bifurcation, extending understanding of oscillation emergence in complex networks.

Findings

Feedback cycles of length three or more support oscillations.

Examples include chemical, ecological, metabolic, and genetic networks.

Reaction kinetics are not limited to standard models.

Abstract

Autonomous sustained oscillations are ubiquitous in living and nonliving systems. As open systems, far from thermodynamic equilibrium, they defy entropic laws which mandate convergence to stationarity. We present structural conditions on network cycles which support global Hopf bifurcation, i.e. global bifurcation of non-stationary time-periodic solutions from stationary solutions. Specifically, we show how monotone feedback cycles of the linearization at stationary solutions give rise to global Hopf bifurcation, for sufficiently dominant coefficients along the cycle. We include four example networks which feature such strong feedback cycles of length three and larger: Oregonator chemical reaction networks, Lotka-Volterra ecological population dynamics, citric acid cycles, and a circadian gene regulatory network in mammals. Reaction kinetics in our approach are not limited to mass action or Michaelis-Menten type.

Equations

$$\dot{x}=f(x)=\sum _j({\bar{y}}_j-y_j)r_j(x)$$

$$r_j>0.$$

$$j\text{mathchar}58y_{j1}{X}_1+\dots +y_{jM}{X}_M⟶{\bar{y}}_{j1}{X}_1+\dots +{\bar{y}}_{jM}{X}_M.$$

$${\bar{y}}_{jm}>y_{jm}.$$

$$r_j(x)=k_j{x}^{y_j}\text{mathchar}58=k_j{x}_1^{y_{j1}}\cdot \dots \cdot x_M^{y_{jM}},$$

$$r_j(x)=k_j\prod _m(x_m/(1+c_{jm}{x}_m){)}^{y_{jm}}$$

$$r_{jm}\text{mathchar}58={\partial }_{x_m}{r}_j>0⟺y_{jm}>0.$$

$$0=f(x^{\ast })=\sum _j({\bar{y}}_j-y_j)r_j(x^{\ast }),$$

$$V(x)\text{mathchar}58=\sum _m{f}_m(x)\cdot v(x_m/x_m^{\ast }),$$

$$\sum _{j\text{mathchar}58{\bar{y}}_j=y}{r}_j(x^{\ast })=\sum _{j^{\prime }\text{mathchar}58y_{j^{\prime }}=y}{r}_{j^{\prime }}(x^{\ast }),$$

$$\prod _{j=1}^N{k}_j^{+}=\prod _{j=1}^N{k}_j^{-},$$

$${\dot{x}}_m=f_m(x_{m-1},x_m,x_{m+1}),$$

$$f_m=k_{m-1}^{+}{x}_{m-1}-(k_{m-1}^{-}+k_m^{+})x_m+k_m^{-}{x}_{m+1}$$

$${\dot{x}}_m=f_m(x),$$

$$f_{m{m}^{\prime }}\text{mathchar}58={\partial }_{x_{m^{\prime }}}{f}_m(x^{\ast }),$$

$${\beta }_k\text{mathchar}58=f_{m_k{m}_{k-1}}\ne 0$$

$$\beta \text{mathchar}58=\prod _{k=1}^N{\beta }_k\ne 0.$$

$$a_k\text{mathchar}58=-f_{m_k{m}_k}\ne 0,$$

$$0=f(\epsilon ,a,x^{\ast })$$

$$f_{x}(\varepsilon,a,x^{*})=\bigg{(}\begin{array}[]{cc}\mathbf{A}+\varepsilon\mathbf{A}^{\prime}&\varepsilon\mathbf{B}\\ \varepsilon\mathbf{C}&\varepsilon\mathbf{D}\end{array}\bigg{)},$$

$$f_{m{m}^{\prime }}=0\text{for}1\le m,m^{\prime }\le N,\text{unless}{m}^{\prime }\in \{m,m-1\}(\mathrm{mod}N).$$

$$\dot{\xi}=\mathbf{A}\xi,\qquad\mathbf{A}=\left(\begin{array}[]{cccc}-a_{1}&&&\beta_{1}\\ \beta_{2}&-a_{2}&&\\ &\ddots&\ddots&\\ &&\beta_{N}&-a_{N}\end{array}\right),$$

$${\beta }_2=\dots ={\beta }_N=+1\text{and}$$

$$\beta =\prod _{m=1}^N{\beta }_m={\beta }_1=\pm 1;$$

$$a_{m}(a)=a\alpha_{m},\qquad\prod_{m=1}^{N}\alpha_{m}=(-1)^{N_{\mathrm{aut}}}\,,\quad a^{N}=\mathrel{\bigg{|}}\prod_{m=n}^{N}a_{m}/\prod_{m=1}^{N}\beta_{m}\mathrel{\bigg{|}}\,,$$

$$⟨\alpha ⟩$$

$$⟨1/\alpha ⟩$$

$$⟨\alpha {⟩}_h$$

$$N_{\mathrm{aut}}-2\lfloor N/4\rfloor -1$$

$$|N_{\mathrm{aut}}-2\lfloor N/4\rfloor -1|$$

$$N_{\mathrm{aut}}-2N/4+\sigma$$

$$|N_{\mathrm{aut}}-2N/4+\sigma |$$

$$N_{\mathrm{aut}}-2\lfloor (N+2)/4\rfloor$$

Full text

We include four example networks which feature such strong feedback cycles of length three and larger: Oregonator chemical reaction networks, Lotka-Volterra ecological population dynamics, citric acid cycles, and a circadian gene regulatory network in mammals. Reaction kinetics in our approach are not limited to mass action or Michaelis-Menten type.

Institut für Mathematik

Freie Universität Berlin

Arnimallee 3

14195 Berlin, Germany

1. 1 Reaction networks and oscillations

2. 2 Main result

3. 3 Global Hopf bifurcation

4. 4 Linear feedback cycles

5. 5 Main result: proof

6. 6 Four examples

7. 6.1 Oregonators

8. 6.2 Lotka-Volterra networks

9. 6.3 Citric acid cycles

10. 6.4 Circadian gene regulation

1 Reaction networks and oscillations

Network graphs are a common modeling device to describe dependencies of certain sub-units among each other. Vertices indicate those sub-units. Directed edges indicate coupling directions, or positive and negative signs of influence. Popular examples in a differential equations context are chemical reaction systems, neural networks, power grids, and many others. Where emphasis may have been on equilibration and steady state behavior, originally, more recent focus has shifted much towards the complexities of temporal and spatial patterns of the collective vertex behavior. The main objective, in the present paper, is to explore the potential of network structures, as such, towards autonomously time periodic network responses. Particulars of coupling parameters will play a subordinate role in that quest. Mostly we address large ranges of parameters. A slow-fast constraint, however, will emphasize a select feedback cycle in the network.

Let us be more specific. Chemical reaction networks, for example, take the form

[TABLE]

with positive vectors $x$ of the metabolite concentrations $x_{m},\ m=1,...,M,$ finitely many nonnegative stoichiometric coefficient vectors $y_{j}\neq\bar{y}_{j}\in\mathbb{R}^{M}$, and positive reaction rate functions

[TABLE]

In chemical notation, the $j$-th summand in (1.1) accounts for the reaction

[TABLE]

One possibility to view (1.1) as a network takes the metabolites $X_{m}$ with concentrations $x_{m}$, as vertices $m$, and dependencies of $f_{m}$ on $x_{m^{\prime}}$, as directed edges $m^{\prime}\rightarrow m$. See our comments on the general setting (2.1) below. Another possibility, suggested by (1.3), is to take the vectors $y,\bar{y}$ as vertices, with reaction arrows as edges. See (1.10) below.

Educts or inputs $m$ of reaction $j$ are defined by nonzero $y_{jm}>0$, and outputs by nonzero $\bar{y}_{jm}$. Nonzero $\bar{y}_{jm}=y_{jm}$ describe a catalyst $m$, for which reaction $j$ does not affect $x_{m}$. Strong autocatalysis of $m$, which catalyzes its own net production, is described by

[TABLE]

Mass action kinetics, prevalent in large parts of classical anorganic chemistry, and in gas phase reactions in particular, is defined by

[TABLE]

usually for integer-valued $y_{j},\bar{y}_{j}$, with the convention $x_{m}^{0}\mathrel{\mathop{\mathchar 58\relax}}=1$. The rate coefficients $k_{j}$ are assumed to be strictly positive; see (1.2). Reactions catalyzed by enzymes, ubiquitous in biological metabolic networks, allow for more general Michaelis-Menten kinetics of the form

[TABLE]

with saturation coefficients $c_{jm}\geqslant 0$. Usually $y_{jm}\in\{0,1\}$. Note how (1.6) reduces to mass action (1.5), for $c_{jm}=0$. The denominator is often “linearized”, inappropriately, to become $1+c^{T}_{j}x$. The partial derivatives of (1.6) satisfy

[TABLE]

Enzymatic inhibition of $r_{j}$ by $x_{m}$, in contrast, is characterized by $\bar{y}_{jm}=y_{jm}$ and a factor $1/(1+c_{jm}x_{m})$, so that $r_{jm}<0$ instead. See sections 6.3 and 6.4 for examples. Therefore we will not subscribe to any monotonicity constraints, in the general mathematical setting of section 2 below. For in-depth information on the very rich subject of chemical reaction kinetics we refer to the currently 43 volumes of the book series Comprehensive Chemical Kinetics [CCK]. For a comprehensive background on chemical reaction networks see [Fei19], by a leading pioneer in the field.

Thermodynamics of closed systems advocates convergence to steady state equilibria $x^{*}$ of (1.1), i.e. to stationary solutions

[TABLE]

The basic tool is a Lyapunov function

[TABLE]

with $v(\xi)\mathrel{\mathop{\mathchar 58\relax}}=\xi\log\xi-\xi+1$, related to the negative of relative entropy. The defining property of Lyapunov functions is that $t\mapsto V(x(t))$ is decreasing along solutions $x(t)$, usually strictly when $x(t)$ is nonstationary. The Lyapunov property of $V$ has been established in [HJ74] under the assumptions of mass action kinetics (1.5) and the following complex balance condition (1.10).

From (1.3) we recall how the stoichiometric vectors $y_{j},\bar{y}_{j}$ may be taken as vertices of the complex graph $\mathcal{C}$, possibly including the complex $y_{j}=0$ and/or $\bar{y}_{j}=0$. Note how identical vectors $y_{j}$ or $\bar{y}_{j}$ for different $j$ may describe the same vertex complex. The directed edges $j$ of $\mathcal{C}$ are simply the reaction arrows (1.3) of standard chemical notation. Then complex balance requires the existence of a positive equilibrium $x^{*}>0$ such that the inflows and outflows balance, i.e.

[TABLE]

at every nonzero complex $y$. In other words, the total production rate at any complex $y=\bar{y}_{j}$ as an output of reactions $j$, balances the total consumption rate at the same complex $y=y_{j^{\prime}}\,$, as an input of other reactions $j^{\prime}\neq j$, just as currents do at vertices of Kirchhoff circuits.

Detailed balance in reversible reaction systems $j^{\pm}$: $y_{j^{\pm}}\leftrightharpoons\bar{y}_{j^{\pm}}$, where $\bar{y}_{j^{+}}=y_{j^{-}}$ and $\bar{y}_{j^{-}}=y_{j^{+}}$ for all $j^{\pm}$, is a special case of complex balance, already considered by [Weg1902]. It requires $r_{j^{+}}(x^{*})=r_{j^{-}}(x^{*})$ for every reversible pair $j^{\pm}$. For reversible monomolecular cycles $y_{j^{+}}=\mathbf{e}_{j},\ y_{j^{-}}=\mathbf{e}_{j+1}\,,\ j\,(\mathrm{mod}\ N)$ and mass action kinetics, detailed balance amounts to the famous Wegscheider condition

[TABLE]

which prevents oscillations. Wegscheider’s arguments for (1.11) were based on thermodynamic considerations on irreversibility, at the microscopic level. It is a lasting merit of [Hir1911] to point at the possibility of (damped) oscillations, once the Wegscheider constraints (1.11) are strongly violated. Our emphasis below on unidirectional $N$-cycles, as a cause for global Hopf bifurcation, is essentially based on this insight.

In passing we note how reversible monomolecular cycles lead to Jacobi systems

[TABLE]

for $m\,(\mathrm{mod}\ N)$ with strictly positive off-diagonal partial derivatives $\partial_{x_{m\pm 1}}f_{m}$. See [FuOl88] for a detailed study. Standard mass action makes

[TABLE]

linear, and $x_{1}+\ldots+x_{N}\equiv\mathrm{const}$ is preserved. Spectral analysis, similar to the case $\beta=+1$ in proposition 4.1 below, then implies stability of steady states, due to the presence of a positive (left) kernel vector. Alternatively, complex balance for the unit vector complexes $y_{m}=\mathbf{e}_{m}$ can be invoked. Note how the addition of strongly autocatalytic diagonal terms like $X_{m}\rightarrow 2X_{m}$ can lead to sustained oscillations and instability. Similar remarks apply to $N$-cycles with general monotone reaction rates $r^{\pm}_{j}$.

Complex balance is clearly sufficient for $x^{*}$ to be a stable steady state (1.8) of (1.1). Notably [Mie17] has much extended the ODE stability results, for the mass action Lyapunov function $V$ in (1.9), to a reaction-diffusion PDE context under Neumann boundary conditions. His results start from a general observation in [Ali79]. Mielke’s extensions include exponential convergence results and cover the presence of stoichiometric invariant subspaces.

Our present paper will remain in an ODE setting, for simplicity of presentation, even though our approach extends to the PDE setting of reaction-diffusion systems. Based on fast $N$-cycles, we study the appearance of time periodic solutions, in contrast to equilibration and beyond the variational complex balance setting (1.10).

Experimental evidence for chemical oscillations has become overwhelming, by now [Zha91]. We recall a few highlights. As early as 1828, Fechner has observed transient polarity reversals in an electro-chemical experiment [Fe1828]; see also [He1901]. The celebrated integrable Lotka-Volterra model [Lot1920] has been described by Lotka, originally, as a hypothetical model for sustained time-periodicity in a chemical reaction with mass action kinetics, and not in the tradited Volterra context [Vol1931] of predator-prey population dynamics. A first chemical experiment with sustained autonomous oscillations was described in [Bray1921]. Experiments on the now famous Belousov-Zhabotinsky reaction (BZ) by Belousov in the 1950s were rejected, at first, on “obvious” thermodynamic grounds. The BZ reaction is an open system, in fact, and the eventual process of dieing out is irrelevant over the time-span of its oscillations. A decade later, young Zhabotinsky rehabilitated the findings by Belousov, and managed to get published [Zha64]. The famous Brusselator “model” [Lef68, PriLef68] for the BZ reaction, by Prigogine and co-workers, had originally just been designed to exhibit and numerically investigate Turing instability [Tur52]. A model for observed glycolytic oscillations in the metabolism of yeast cells was suggested by [Sel68]. All the above considerations were based on phase plane analysis, i.e. on reaction systems (1.1) with $M=2$ metabolites $m=1,2$. The article [Hig67] provides a comprehensive survey and discussion of the planar possibilities.

The chemically more realistic Oregonator model [FN74] of the BZ-reaction was an early example of oscillatory chemical reaction networks involving at least $M=3$ metabolites; see section 6.1. Eigen’s quite hypothetical hypercycle [Eig71], of course, also known as the replicator equation, features cycles of any length $N$ in an attempt to model molecular evolution; see also the book [HS98]. It can be seen as a projective version, for population percentages, of general Lotka-Volterra models [Oli14] discussed in section 6.2. Oscillations in the famous citric acid cycle (CAC, Krebs cycle) involving eight metabolites have been described, experimentally, by [MacDetal03]; see section 6.3 below. In section 6.4 we discuss a gene regulatory network for circadian rhythms in mammals [Miretal09]. Non-isothermal oscillations, where the temperature dependence of the rate functions $r_{j}$ plays a decisive role, have been studied much, in the PDE context of spatially heterogeneous catalysis. See [Aris75, Fie83] for experimental and mathematical results, as well as [IE95] for a survey of the early developments.

Theoretical results on autonomous time-periodic oscillations are rare, beyond mere numerical simulation. Mostly, they establish the existence of stationary solutions $f(x^{*})=0$ with purely imaginary eigenvalues, by the $M$-dimensional Routh-Hurwitz criterion. Classical local Hopf bifurcation [Hopf1942, MaMcC76, CR77] then is supposed to infer periodic solutions. For mass action kinetics (1.5), however, the computational difficulties seem to grow prohibitively with dimension. Even best analytic results like [GES05, EEetal15], which require advanced techniques and concepts from computational algebra, do not proceed beyond $M=3,4$. They also fail to address standard prerequisites of local Hopf bifurcation, like spectral nonresonance and transverse crossing conditions. Instead, our approach will avoid the restrictions of mass action kinetics, and will explore fast feedback cycles in networks as a source of global Hopf bifurcation.

In the next section we outline our main result on global Hopf bifurcation of periodic orbits induced by fast feedback cycles, in a much more general mathematical setting. See our main result, theorem 2.2. Our precise notion of global Hopf bifurcation involves some subtleties which we postpone to section 3; see definition 3.2 and corollary 3.4. We recall some tools for global Hopf bifurcation there, as developed in [Fie85]; see theorem 3.3 and corollary 3.4. In section 4 we collect the prerequisite spectral properties of cyclic monotone feedback systems, in the spirit of [M-PS90]. In section 5, this enables us to prove theorem 2.2 as an application of corollary 3.4. We conclude with the promised four applications, in section 6.

Acknowledgment. The present results are based entirely on many pleasantly challenging and always inspiring discussions with Alexander Mielke, including his beautiful contribution [Mie17]. The present paper is dedicated to him in deep gratitude and lasting friendship. We are also indebted to Marty Feinberg, for sharing his pioneering insights so generously and encouragingly, over so many years, and to Nicola Vassena for many helpful comments concerning network dynamics. Patricia Hăbăşescu performed the typesetting with outstanding dedication and diligence, and some expert help from Alejandro Lopez. This work was partially supported by the Deutsche Forschungsgemeinschaft through SFB 910 project A4.

2 Main result

In this section we formulate our main result, theorem 2.2, for general ODEs

[TABLE]

$m=1,\ldots,M,$ with $C^{1}$-nonlinearities $f_{m}$ and for $x\in\mathbb{R}^{M}$. Let $x^{*}$ be a stationary solution of (2.1), i.e. $f(x^{*})=0$. Then the Jacobian $f_{x}(x^{*})=(f_{mm^{\prime}})$ is given by the partial derivatives

[TABLE]

for $1\leqslant m,m^{\prime}\leqslant M$.

In network language, the ODE setting (2.1) generalizes (1.1) as follows. Suppose $f_{m}=f_{m}(x_{I(m)})$ . This notation indicates that $f_{m}$ depends on the component $x_{m^{\prime}}$ of $x\in\mathbb{R}^{M}$ if, and only if, $m^{\prime}$ is in the set $I(m)\subseteq\{1,\ldots,M\}$ of inputs of $f_{m}$. Then the dependencies $I(m)$ of $f_{m}$ define a directed graph $\Gamma$ with metabolite vertices $m$ and directed edges $m^{\prime}\rightarrow m$ from $m^{\prime}\in I(m)$ to $m$. We explicitly allow, but do not impose, self-loops $m\in I(m)$. The graph $\Gamma$ in the setting (2.1) is commonly employed in the description of gene regulatory networks; see for example [FieMKS13]. Note $f_{mm^{\prime}}=0$, unless $m^{\prime}\in I(m)$. The full ODE (2.1) corresponds to the maximal bi-directional graph $\Gamma$ with self-loops, of course.

Based on nonzero entries $f_{mm^{\prime}}\,,\ m^{\prime}\in I(m)$, of the Jacobian $f_{x}(x^{*})$, we can now identify fast feedback cycles in the directed graph $\Gamma$, which play the central role for our results on fast oscillations. For the sake of generality, however, we do not restrict ourselves to any specific directed graph $\Gamma$. Rather, we focus on dominant $N$-cycles, just on the linear level of the Jacobian $f_{x}(x^{*})$, as a source of oscillations and of global Hopf bifurcation.

2.1 Definition**.**

Fix $N\in\{2,\ldots,M\}$. Let $\mathbf{m}=(m_{1}\ \ldots\ m_{N})$ denote any ordered $N$-tuple of distinct metabolites $m_{k}\in\{1,\ldots,M\}$. We call $\mathbf{m}$ an $N$-cycle if

[TABLE]

holds, for all indices $k\,(\mathrm{mod}\ N)$. We say the $N$-cycle possesses positive or negative feedback, depending on the sign of

[TABLE]

For a nondegenerate $N$-cycle we require, in addition, nonzero self-loops

[TABLE]

for all $k=1,\ldots,N$. Motivated by reaction network dynamics, we call the number $0\leq N_{\mathrm{aut}}\leq N$ of $a_{k}<0$, i.e. the number of strictly positive self-feedbacks $f_{m_{k}m_{k}}$, the autocatalytic number of the nondegenerate $N$-cycle $\mathbf{m}$.

We can now describe the detailed setting of our main result, theorem 2.2 below. It is of crucial importance here, and deviates significantly from previous work in the area, that we consider the partials $f_{mm^{\prime}}$ as free parameters which may vary independently of the steady state $x^{*}$ and, to some extent, independently of each other. More precisely we consider networks $\dot{x}=f(\varepsilon,a,x)$ depending on a parameter $a>0$ and a small parameter $\varepsilon>0$, such that

[TABLE]

possesses a parameter-independent stationary solution $x^{*}$.

For the Jacobian at $x^{*}$ we assume an expansion

[TABLE]

in block matrix form, with small $\varepsilon>0$. Only $\mathbf{A}=\mathbf{A}(a)$ is allowed to depend on the parameter $a$. Specifically, we assume that the $N\times N$ block matrix $\mathbf{A}=(f_{mm^{\prime}})_{1\leq m,m^{\prime}\leq N}$ of $f_{x}$, at $\varepsilon=0$, describes a nondegenerate $N$-cycle, $\det\mathbf{A}\neq 0$. Here we have relabeled $\mathbf{m}$, without loss of generality, such that $\mathbf{m}=(1\ \ldots\ N)$ in definition 2.1. We require all other entries of $\mathbf{A}$, not supported on the $N$-cycle, to be zero:

[TABLE]

The blocks $\mathbf{A}^{\prime},\mathbf{B},\mathbf{C},\mathbf{D}$ in the $M\times M$ Jacobian (2.7) are assumed to remain bounded, uniformly in the small scaling parameter $\varepsilon>0$. These assumptions make the $N$-cycle described by $\mathbf{A}$ dominant and fast compared to the remaining dynamics, on the linear level.

The mathematical motivation for our emphasis on cycles, in addition to [Hir1911], comes from the Quirk-Ruppert-Maybee theorem; see the beautiful account in [JKD77]. That theorem addresses matrices $\mathbf{M}$ with prescribed sign structure of their entries. It characterizes spectral stability Re spec $\mathbf{M}\leq 0$, for all $\mathbf{M}$ with the same sign structure, by three simultaneous requirements: nonpositive diagonal elements, nonpositive products over 2 cycles, and vanishing products $\beta$ over $N$-cycles, for $N\geq 3$. Our results can be seen as an attempt to assert global Hopf bifurcation when the second conditions is violated, by an autocatalysis count $N_{\mathrm{aut}}>0$, and/or the third condition is violated, by fast $N$-cycles $\mathbf{A}$ with $\mathrm{sign}\ \beta=\pm 1$.

The linearization

[TABLE]

on the fast $N$-cycle constitutes a linear cyclic monotone feedback system. See [M-PS90] for a detailed spectral analysis and deep nonlinear consequences. For simplicity we perform linear rescalings of each $\xi_{m}$ , and a positive linear rescaling of time $t$, to normalize the off-diagonal $N$-cycle elements of $\mathbf{A}$ such that

[TABLE]

see (2.4). Note here that the product $\beta$ of the $\beta_{m}$ is invariant under scalings of the $\xi_{m}$ . Under positive time rescaling the nonzero feedback sign becomes $\mathrm{sign}\,\beta=\beta=\beta_{1}=\pm 1$. We define the bifurcation parameter $a>0$ by the normalization

[TABLE]

along the original diagonal of the nondegenerate $N$-cycle $\mathbf{A}$; see (2.5). For the normalized diagonal elements $\alpha_{m}\neq 0$ we use the abbreviations $\langle\cdot\rangle$ and $\langle\cdot\rangle_{h}$ to denote their arithmetic and harmonic means, respectively. We assume

[TABLE]

whenever $\sigma,\sigma_{h}$ appear.

In addition to the signs of the above arithmetic and harmonic means, our oscillation conditions will only involve the length $N\geqslant 3$ of the catalytic cycle, and the count $N_{\mathrm{aut}}$ of diagonal strongly autocatalytic entries $\alpha_{m}<0$. Specifically, we assume any one of the following four cases to hold.

1. (i)

For positive feedback $\beta=+1$ and $N\not\equiv 0\ (\mathrm{mod}\ 4)$:

[TABLE] 2. (ii)

For positive feedback $\beta=+1$ and $N\equiv 0\ (\mathrm{mod}\ 4)$:

[TABLE] 3. (iii)

For negative feedback $\beta=-1$ and $N\not\equiv 2\ (\mathrm{mod}\ 4)$:

[TABLE] 4. (iv)

For negative feedback $\beta=-1$ and $N\equiv 2\ (\mathrm{mod}\ 4)$:

[TABLE]

Here $\lfloor\cdot\rfloor$ denotes the integer floor function. We repeat that the signs $\sigma,\sigma_{h}=\pm 1$ of the arithmetic and harmonic means are assumed to be nonzero, respectively, whenever they appear. Specifically, this assumption will only arise under the following circumstances:

[TABLE]

Concerning the open parameter interval $a\in(\underline{a},\bar{a})$ where we will assert global Hopf bifurcation, we distinguish the following cases. We fix

[TABLE]

with arbitrarily small $\underline{a}$ in the second case. Similarly, we fix

[TABLE]

with arbitrarily small $1-\bar{a}$ in the second case.

2.2 Theorem**.**

Consider a network (2.1) with a fast nondegenerate $N$-cycle $\mathbf{A}$ on $\mathbf{m}=(1\ \ldots\ N),\ N\geq 3$, of the Jacobian (2.7) – (2.9) at the $(\varepsilon,a)$-independent steady state $x^{*}$ in (2.6). Assume hyperbolicity of the lower right block $\mathbf{D}$ in the Jacobian (2.7), i.e. $0\,\not\in\,\mathrm{Re}\,\mathrm{spec}\,\mathbf{D}$. Let the normalization (2.10) and parameter assumptions (2.11), (2.17) – (2.19), on $a>0$ hold.

Then each of the cases (2.13) – (2.16) leads to the following conclusion. There exists $\varepsilon_{0}>0$ depending on $\underline{a},\bar{a}$ such that for any fixed $0<\varepsilon<\varepsilon_{0}$ the network (2.1) exhibits global Hopf bifurcation of nonstationary periodic solutions from the steady state $x^{*}$, for parameters $a\in(\underline{a},\bar{a})$.

For our precise notion of global Hopf bifurcation we refer to definition 3.2 and corollary 3.4 below.

Let us rewrite the four options (2.13) – (2.16) in more concise form. Specifically, we define the abbreviating signs

[TABLE]

Note the isomorphism $\iota$ from the group $\mathbb{Z}_{4}^{*}$ of multiplicative units in $\mathbb{Z}_{4}$ to the multiplicative group $\mathbb{Z}_{2}=\{\pm 1\}$. The four options (2.13) – (2.16) can then be condensed into the single assumption

[TABLE]

For a first consistency check of our result, consider time reversal $t\mapsto-t$ in our original ODE setting (3.1). Stationary solutions, periodic orbits (as sets), and global Hopf bifurcation remain unaffected by time reversal. But how about the options (2.13) – (2.16) of theorem 2.2, alias assumption (2.22)?

Under time reversal, all $\alpha_{m}\,,\,\beta_{m}$ reverse sign, and so do the signs $\sigma,\sigma_{h}$ of the arithmetic and harmonic means $\langle\alpha\rangle,\langle\alpha\rangle_{h}$ . Sign reversal of all $\alpha_{m}$ therefore replaces the count $N_{\mathrm{aut}}$ of strongly autocatalytic $\alpha_{m}<0$ in the $N$-cycle by $N-N_{\mathrm{aut}}$ . This reverses the sign of the term $2N_{\mathrm{aut}}-N$ in assumption (2.22). The normalized feedback coefficient $\beta$ gets multiplied by $(-1)^{N}$. By abbreviation (2.20), however, which eliminates even $N$, this also just reverses the sign of $\iota\beta$ in assumption (2.22). The coefficient $\kappa=\kappa(N,\beta)$, in contrast, vanishes for odd $N$ and thus remains unaffected by time reversal. In summary, these elementary considerations show how assumption (2.22) remains invariant under time reversal – just as the conclusion of global Hopf bifurcation in theorem 2.2 does.

3 Global Hopf bifurcation

We introduce the main tools in our analysis of autonomous time periodic oscillations. Skipping proofs, we adapt results on global Hopf bifurcation going back to [Fie85], and based on earlier results by Jim Yorke and others [AY78, M-PY82, AllY84]. Specifically, we introduce virtual periods, and the center index

 

$\Box$ (pronounced “zhong”). Our main abstract results are summarized in theorem 3.3 and corollary 3.4 below. See also [Fie88], section 3, for a more detailed survey.

In this section we consider general vector fields

[TABLE]

on $x\in\mathbb{R}^{M}$, with scalar parameter $a\in\mathbb{R}$, and continuous $f,f_{x}$. As before, we call $x^{*}$ stationary at $a$, or steady state, if $f(a,x^{*})=0$. We call a solution $x(t)$ periodic with a period $T>0$ if $x(t)$ is nonstationary and

[TABLE]

holds for all real $t$. The set of all periods $T$ then takes the form $T=kp$, with $k\in\mathbb{N}$, where $p>0$ is called the minimal period of $x(t)$.

Global bifurcation results are often based on topological tools which, in turn, are established via approximations and limit arguments. For periodic solutions, a peculiar difficulty arises from the notion of minimal period $p$, versus periods $T=kp$: the limits of minimal periods may fail to be minimal periods of the limit. At Hopf bifurcation, for example, we observe sequences of periodic orbits shrinking to a stationary solution, although the minimal periods approach a positive limit. At period doubling, we observe sequences of periodic orbits where the minimal period drops by a factor two, in the limit. More generally, the same difficulty arises in systems with group equivariance: isotropies may increase (but not decrease) under limits. This more encompassing equivariant viewpoint, which we do not pursue any further here, has been addressed in [Fie88], section 4.

The notion of virtual periods remedies the limit deficiency of minimal periods, by including the associated linear flow. We call $q>0$ a virtual period of $x(t)$ at $a$, if $q$ is the minimal period of some pair $(x(t),y(t))$, where $y(t)$ satisfies the (nonautonomous) linearized equation

[TABLE]

for all real $t$. We also use this terminology if $x(t)\equiv x^{*}$ happens to be stationary. In summary, virtual periods are the minimal periods of the induced flow on the tangent bundle.

Stationary $x^{*}$ with virtual periods $q$ are called Hopf points: indeed they possess nonzero purely imaginary eigenvalues. Standard Hopf bifurcation, from a transverse crossing of a pair of simple and nonresonant eigenvalues $\pm\mathrm{i}$ of $f_{x}(a,x^{*})$, is indicative of a virtual period $q=2\pi$ at $x^{*}$. Any subset of rationally related eigenvalues $\pm\mathrm{i}\omega$ generates a virtual period $q$ given by the least common multiple of the $2\pi/\omega$. Standard period doubling at minimal period $p$ features two virtual periods $q$: $p$ itself, and $2p$. Note how the set of virtual periods $q$ of any fixed stationary or periodic solution $x(t)$ is always bounded, whereas the set $T=kp$ of all periods is always unbounded.

It turns out that virtual periods, unlike minimal periods, are closed under limits.

3.1 Proposition**.**

Let $q_{n}$ be a virtual period of $x_{n}$ at parameter $a_{n}$. Assume bounded convergence:

[TABLE]

Then $q_{\infty}>0$, and $q_{\infty}$ is a virtual period of $x_{\infty}$ at parameter $a_{\infty}$.

Henceforth we require all Hopf points $(a_{*},x^{*})$ of (3.1) to be nondegenerate, i.e.

[TABLE]

This allows us to continue the steady state $x^{*}=x^{*}(a)$, locally, by the implicit function theorem. Let

[TABLE]

count the strictly unstable eigenvalues at $(a,x^{*}(a))$, with algebraic multiplicity. We now require Hopf points to be isolated, in $\mathbb{R}\times\mathbb{R}^{M}$. Then $\mu(a)$ is the unstable dimension, or Morse index, of the hyperbolic steady state $x^{*}(a)$, for nearby a\neq a_{*}\. This allows us to define the crossing number

[TABLE]

of the Hopf point $x^{*}$ at $a=a_{*}$ . This is the net number of eigenvalue pairs crossing the imaginary axis from left to right, as $a$ increases through $a_{*}$ . Finally, following [M-PY82], we define the center index of the Hopf point $(a_{*},x^{*})$ as

[TABLE]

Fix any open subset $\mathcal{U}$ of $\mathbb{R}\times\mathbb{R}^{M}$, such that $\mathcal{U}$ contains the whole nonstationary periodic orbit, with any point on it. We clarify our notion of global versus local continua of periodic solutions and Hopf points in $\mathcal{U}$ as follows. Denote

[TABLE]

In other words, $\mathcal{P}=\check{q}\tilde{Q}$, where the projection $\check{q}$ omits the $q$-component of $\mathcal{Q}$.

3.2 Definition**.**

A connected component $\mathcal{C}$ of $\mathcal{P}$, i.e. of the periodic solutions and Hopf points in $\mathcal{U}$, is called local in $\mathcal{U}$, if the closure of $\mathcal{C}$ is compactly contained in $\mathcal{U}$ and the virtual periods in $\mathcal{C}$ are bounded above. In other words, the lift $\check{q}^{-1}\mathcal{C}$ is compactly contained in $\mathcal{Q}$. Connected components $\mathcal{C}$ which are not local are called global.

Note that proposition 3.1 asserts compactness of local components $\mathcal{C}$.

3.3 Theorem**.**

Consider the flow (3.1) and assume all Hopf points in $\mathcal{U}$ are nondegenerate, as in (3.5), and isolated. Let $\mathcal{C}$ be a local connected component of the periodic solutions and Hopf points $\mathcal{P}$ in $\mathcal{U}$. Then

[TABLE]

where the sum ranges over the finitely many Hopf points in $\mathcal{C}$, if any.

3.4 Corollary**.**

In the setting of theorem 3.3, assume $\mathcal{P}$ contains only finitely many Hopf points and

[TABLE]

Then $\mathcal{P}$ possesses at least one global connected component $\mathcal{C}$ which also contains a Hopf point. We call this case global Hopf bifurcation in $\mathcal{U}$.

To derive the corollary from the theorem, let $\mathcal{C_{\ell}}$ enumerate the finitely many disjoint connected components of $\mathcal{P}$, which contain Hopf points. Suppose, indirectly, that each $\mathcal{C_{\ell}}$ is local. Then (3.10) implies

[TABLE]

contradicting (3.11). Hence at least one $\mathcal{C}_{\ell}$ is global, by theorem 3.3.

The proof of theorem 3.3 is based on generic approximation. The cancellation (3.10) of center indices on compact connected components $\mathcal{C}$ follows from the same property in the generic situation, by approximation. See [M-PY82] for the generic case. This requires a parametrized version of the Kupka-Smale theorem, via Thom-transversality, and a detailed degree argument which carefully distinguishes periodic orbits with orientable and nonorientable unstable manifolds. The resulting global Hopf components of orientable periodic orbits in $\mathcal{P}$ were called snakes, in [M-PY82]. The only discontinuities of minimal periods, in the generic case of snakes, occur at period doubling bifurcations. By generic approximation, this reveals jumps by factors 2 as the only possible discontinuities of virtual periods, in the general case of nongeneric snake limits.

Competing topological results, based on the $J$-homomorphism or $S^{1}$-equivariant degree theory, studied continua of triples $(a,x,T)$ with (not necessarily minimal) periods $T$ of $(a,x)$ as in (3.2); see [AY78] for the original, and [IV03] for more recent developments with many references. The intriguing “jug-handle” example by [AllY84] exhibits a continuum with bounded $(a,x)$ and unbounded $T$, whereas virtual periods remain bounded. The jug-handle consists of a compact loop of periodic orbits $(a,x)$, where the two branches generated at a saddle-node bifurcation re-unite, at a period-doubling. Any such loop generates an unbounded continuum of triples $(a,x,T)$ where $T$ traverses all multiples $2^{k}p$ of the minimal periods $p$, for $k\in\mathbb{N}_{0}$. The jug handle, which is evidently quite bounded in parameter and phase space, therefore counts as a global, unbounded continuum in the sense of [AY78, IV03]. The virtual periods on the jug handle, however, remain bounded: they are given by $p$ and, at the period doubling only, by $\{p,2p\}$. The same jug handle therefore does not qualify as global, in the sense of [M-PY82, Fie85] which we adopt in the present paper.

Of course our notion of “globality” depends on the choice of the underlying open domain $\mathcal{U}\subseteq\mathbb{R}\times\mathbb{R}^{M}$ where we study our continua. Indeed we can only assert a global trichotomy for any global component $\mathcal{C}\subseteq\mathcal{U}$:

(3.13)

* *

(i)

either $\mathcal{C}$ is unbounded, or

(ii)

$\mathcal{C}$ is bounded, but $\partial\mathcal{C}\cap\partial\mathcal{U}\neq\emptyset$, or else

(iii)

$\mathrm{clos}\ \mathcal{C}$ is compactly contained in $\mathcal{U}$ with unbounded virtual periods.

Option (iii) of the global trichotomy (3.13) is particularly interesting. For example, consider a convergent sequence $(a_{n},x_{n})\rightarrow(a_{\infty},x_{\infty})$ of Hopf points with purely imaginary eigenvalues $\pm i\omega_{n}$, such that $\omega_{n}\searrow 0$. The steady state $x_{\infty}$ then features an eigenvalue $\omega_{\infty}=0$ with algebraic multiplicity at least two; in the simplest interesting case this is a Bogdanov-Takens point. The virtual periods $q_{n}\mathrel{\mathop{\mathchar 58\relax}}=2\pi/\omega_{n}$ converge to $+\infty$, of course.

More generally, suppose $(a_{n},x_{n})$ are nonstationary periodic with minimal periods $p_{n}\rightarrow+\infty$. Suppose $(a_{n},x_{n})\rightarrow(a_{\infty},x_{\infty})$ becomes stationary, in the limit, but some part of the periodic orbit $x_{n}(t)$ of $x_{n}$ does not converge to $x_{\infty}$. In the simplest interesting case this may happen by convergence of $x_{n}(.)$ to a homoclinic orbit attached to the steady state $x_{\infty}\,$. This example is closely related to the Takens-Bogdanov case, which generates small amplitude homoclinic orbits. For global consequences in vector fields with two real parameters see [Fie86, Fie96].

Suppose next that the orbits $x_{n}(.)$ remain bounded, and stay away from any steady states. Remarkably $p_{n}\rightarrow\infty$ can still occur, along a continuum of periodic orbits and without any bifurcations affecting the minimal periods $p_{n}$. In 1995, such blue sky catastrophes have first been constructed by Turaev and Shilnikov, in a structurally stable way involving a single parameter. See the survey [SST14].

In section 4, we will apply corollary 3.4 to the situation of theorem 2.2. In particular we note how the crossing numbers in (3.11) simply add up to a net crossing number, along a steady state $x^{*}$ which does not actually depend on $a$, as long as $x^{*}$ remains nondegenerate. To account for the slow-fast dichotomy (2.7) of the linearization $f_{x}(\varepsilon,a,x^{*})$ we will also narrow attention from $a\in(-\infty,+\infty)$ to $0<a\in(\underline{a},\bar{a})$.

4 Linear feedback cycles

In this section we collect some spectral properties of the normalized nondegenerate $N$-cycle

[TABLE]

with $a>0,\ \prod\alpha_{m}=(-1)^{N_{\mathrm{aut}}}\,$, and $\beta=\pm 1$. See (2.7) – (2.12) and theorem 2.2.

Proposition 4.1 recalls the general pairwise ordering of the eigenvalues $\lambda_{k}$ of $\mathbf{A}$ by their real parts, due to [M-PS90]. Proposition 4.2 addresses crossings of eigenvalues through the imaginary axis, at $a=0$ and $a=1$. Proposition 4.3 collects the limits, at $a=0^{+},1^{\pm}$, and $\infty$, of the strict unstable dimensions $\mu(a)$ introduced in (2.12):

[TABLE]

We conclude in proposition 4.4, by showing how the presence of a zero eigenvalue, at $a=1$, prevents all further purely imaginary eigenvalues to occur for any $a\geq 1$.

The zero number $z(\xi)$, an integer-valued Lyapunov function for $\dot{\xi}=\mathbf{A}\xi$, is the crucial tool in the deep spectral (and nonlinear) analysis of [M-PS90]. In our normalization (4.1), consider positive feedback $\beta=+1$ first and let $0\neq\xi\in\mathbb{R}^{N}$. Then $z(\xi)\geq 0$ denotes the (even) number of strict sign changes in the ordered cyclic sequence of $\xi$-components $\xi_{m}$, with $m\,(\mathrm{mod}\ N)$. For negative feedback $\beta=-1$, however, we modify that count between $\xi_{N}$ and $\xi_{1}$, only, to account for a strict sign change between $\beta\xi_{N}=-\xi_{N}$ and $\xi_{1}$, instead. In particular $z(\xi)\geq 1$ becomes odd. In summary we obtain

[TABLE]

for both feedback cases, $\beta=\pm 1$.

4.1 Proposition**.**

Assume negative feedback, $\beta=-1$. Then the eigenvalues $\lambda_{k}$ of $\mathbf{A}$ can be ordered in pairs, repeated with algebraic multiplicity, such that

[TABLE]

for indices ranging from [math] to $N-1$. Simple real eigenvalues are labeled in strictly decreasing order. In particular, all eigenvalues are algebraically simple, except for some possibly double real eigenvalues.

The associated real eigenvectors $\xi_{k}$ of $\lambda_{k}$ can be chosen to satisfy

[TABLE]

Here $\xi_{2k}$ and $\xi_{2k+1}$ refer to the real and imaginary parts of the complex eigenvectors, in case $\lambda_{2k+1}=\bar{\lambda}_{2k}$ are nonreal conjugate complex.

For positive feedback $\beta=+1$, the analogous ordering reads

[TABLE]

with real eigenvectors $\xi_{k}$ of $\lambda_{k}$ satisfying

[TABLE]

for resulting indices in $\{0,\ldots,N-1\}$.

Proof..

See [M-PS90]. ∎

To get slightly more specific we write the characteristic equation for the characteristic polynomial $\mathbf{p}$ of $\mathbf{A}$ from (4.1) as

[TABLE]

Here we have used the normalizations (2.11), and the notation (2.12) for arithmetic means. The case $a=0$ in (4.8), with the $N$-th roots of unity $\lambda_{k}^{N}=\beta=\pm 1,\ k=0,\ldots,N-1$, as simple eigenvalues, provides an instructive example for the two feedback cases of proposition 4.1.

4.2 Proposition**.**

Consider the normalized nondegenerate $N$-cycle $\mathbf{A}(a)$ of (4.1), for $a\geq 0$. Then the following holds true.

1. (i)

An eigenvalue $\lambda_{k}=0$ occurs if, and only if,

[TABLE]

The eigenvector $\xi_{k}$ of $\lambda_{k}=0$ satisfies

[TABLE]

The eigenvalue $\lambda_{k}=0$ is simple if, and only if,

[TABLE]

In that case, $\lambda_{k}(a)$ crosses the imaginary axis transversely, at $a=1$, with nonzero derivative

[TABLE]

given by the harmonic mean; see (2.12). 2. (ii)

At $a=0$, the eigenvalues $\lambda_{k}(a)$ are given by the $N$ simple roots of unity

[TABLE]

with $k=0,\ldots,N-1$. Their derivatives with respect to $a$, at $a=0$, are all equal, given by the arithmetic mean

[TABLE]

In particular, the purely imaginary eigenvalues $\lambda_{k}=\pm i$ which occur for $N\equiv 0,2\,(\mathrm{mod}\ 4)$ and $\beta=+1,-1$, respectively, cross the imaginary axis transversely, for arithmetic means $\langle\alpha\rangle\neq 0$.

Proof..

We use expansion (4.8) of the characteristic polynomial.

To prove claim (4.9) of (i), we just insert $\lambda=0$ in (4.8) and recall $a\geq 0$. Algebraic simplicity claim (4.11) is equally obvious, for $\langle 1/\alpha\rangle\neq 0$. Implicit differentiation of $\mathbf{p}(a,\lambda(a))=0$ at $a=1,\ \lambda(1)=0$ yields

[TABLE]

which proves (4.12).

To prove (4.10) note that any eigenvector $0\neq\xi\in\ker\mathbf{A}$ of $\lambda=0$ at $a=1$ satisfies

[TABLE]

for $m\,(\mathrm{mod}\ N),\ \beta_{1}=\beta$, and $\beta_{2}=\ldots=\beta_{N}=1$; see (4.1). Consider the case $\beta=+1,\ \alpha_{1}>0$, first. Then those $2\leq k\leq N$ with strongly autocatalytic $\alpha_{k}<0$ indeed provide precisely $N_{\mathrm{aut}}$ sign changes in the cyclic sequence of $\xi_{m}$, for $m\,(\mathrm{mod}\ N)$. The remaining cases of (4.10) are proved analogously. This proves claim (i).

To prove claim (ii), we insert $a=0$ in $\eqref{eq:3.8}$ to obtain the algebraically simple $N$-th roots of unity $\lambda_{k}^{N}=\beta,\ k=0,\ldots,N-1$. Implicit differentiation of $\mathbf{p}(a,\lambda(a))=0$ at $a=0,\ \lambda^{N}=\beta\neq 0$, indeed yields

[TABLE]

This completes the proof of the proposition. ∎

4.3 Proposition**.**

At $a=\infty$ we obtain the following limiting strict unstable dimension:

[TABLE]

Assume any $\mathbf{A}(a),\ a\geq 0$, possesses an eigenvalue $\lambda_{k}=0$, in the ordering of proposition 4.1. Then $(-1)^{N_{\mathrm{aut}}}=\beta$ and $a=1$, by (4.9).

If we also assume $\langle 1/\alpha\rangle\neq 0$, as in (2.12), so that the harmonic mean with sign $\sigma_{h}=\mathrm{sign}\langle\alpha\rangle_{h}=\pm 1$ exists, then the limiting strict unstable dimensions $\mu(1^{\pm})$ in (4.2) are

[TABLE]

and the even/odd parity of $k$ is determined by

[TABLE]

In the limit $a\searrow 0$ and for $N\not\equiv 1-\beta\ (\mathrm{mod}\ 4),\ \beta=\pm 1$, we obtain

[TABLE]

For $N\equiv 1-\beta\ (\mathrm{mod}\ 4)$ we assume $\sigma=\mathrm{sign}\langle\alpha\rangle\neq 0$, as in (2.12), and obtain

[TABLE]

Proof..

To prove $\mu(\infty)=N_{\textrm{aut}}$ we invoke the characteristic equation (4.8), once again. Trivially, $-\alpha_{m}\neq 0$ with $m=1,\ldots,N$ enumerate the limits of $\lambda_{k}(a)/a$ with $k=0,\ldots,N-1$, for $a\rightarrow+\infty$. This proves claim (4.18).

To prove claim (4.19) we consider the simple eigenvalue $\lambda_{k}=0$ at $a=1$ with eigenvector $\xi_{k}$ and $z(\xi_{k})=N_{\mathrm{aut}}$, from proposition 4.2. In particular, the ordering of $\mathrm{Re}\,\lambda_{k}$ in proposition 4.1 implies $\mu(a)=k$ for the strict unstable dimension $\mu$ at $a=1$. Our assumption $\langle 1/\alpha\rangle\neq 0$ in (2.12) also implies transverse crossing (4.12) of $\lambda_{k}(a)$, at $a=1$, so that

[TABLE]

for small $|1-a|>0$. Because proposition 4.1 excludes any other purely imaginary eigenvalues at $a=1$, besides the simple eigenvalue $\lambda_{k}(a)=0$, this proves claim (4.19).

To prove (4.20) we observe that the absence of zero eigenvalues for $1<a<\infty$ implies that $\mu(1^{+})$ and $\mu(\infty)=N_{\mathrm{aut}}$ share the same parity. Therefore (4.19) implies (4.20) via

[TABLE]

It remains to consider $\mu(0^{+})$ with eigenvalues $\lambda$ at $a=0$ given by the simple roots of unity $\lambda^{N}=\beta=\pm 1$. For the strict unstable dimension $\mu(0)$, which ignores purely imaginary eigenvalues, elementary counting shows $\mu(0)=2\lfloor(N-2+\beta)/4\rfloor+1+(1-\beta)/2$. For $N\not\equiv 1-\beta\ (\mathrm{mod}\ 4)$, purely imaginary roots $\lambda_{k}=\pm i$ do not occur. Therefore $\lfloor(N-2+\beta)/4\rfloor=\lfloor(N-1+\beta)/4\rfloor$ proves (4.21). For $N\equiv 1-\beta\ (\mathrm{mod}\ 4)$, the purely imaginary pair (as all other roots) satisfies $\lambda^{\prime}_{k}=\mathrm{Re}\,\lambda^{\prime}_{k}=-\langle\alpha\rangle$; see (4.14). Therefore

[TABLE]

for small $a>0$, and hence $\mu(0^{+})=\mu(0)+1-\sigma$. Insertion of our elementary count for $\mu(0)$ proves (4.22), and the proposition. ∎

4.4 Proposition**.**

As in (4.19), (4.20) suppose any $\mathbf{A}(a),\ a\geq 0$ possesses an eigenvalue $\lambda_{k}=0$, i.e. $(-1)^{N_{\mathrm{aut}}}=\beta$ and $a=1$. Assume $\langle 1/\alpha\rangle\neq 0$.

Then $\mathbf{A}(a)$ does not possess any other zero or purely imaginary eigenvalues, for any $1\leq a<\infty$, except that simple zero eigenvalue $\lambda_{k}=0$ at $a=1$.

Proof..

Propositions 4.1 and 4.2(i) establish the claim at $a=1$. We have to show that purely imaginary nonzero eigenvalues cannot occur, for any $a>1$. We only consider the case of positive feedback, $\beta=+1$; the case $\beta=-1$ is analogous.

Since $(-1)^{N_{\mathrm{aut}}}=\beta$, positive feedback $\beta=+1$ implies $N_{\mathrm{aut}}$ is even. Therefore (4.10) implies

[TABLE]

for the eigenvector $\xi_{k}$ of $\lambda_{k}=0$ at $a=1$. The ordering (4.6), (4.7) of real eigenvalues implies

[TABLE]

for the two real simple eigenvalues $\lambda_{2k^{\prime}-1}>\lambda_{2k^{\prime}}$, one of them being zero, at $a=1$. We claim $\lambda_{2k^{\prime}-1}$ and $\lambda_{2k^{\prime}}$ straddle zero, for all $a>1$: both eigenvalues remain simple, real, and satisfy

[TABLE]

Then the straddling eigenvalues $\lambda_{2k^{\prime}-1}\,,\lambda_{2k^{\prime}}\,,$ in view of the ordering (4.6), (4.7), prevent any other real or complex eigenvalues from crossing the imaginary axis, at any $a>1$, and the proposition will be proved.

We prove our remaining claim (4.28) for $\sigma_{h}=+1$; the case $\sigma_{h}=-1$ is analogous. Parity property (4.20), $(-1)^{k}=\beta\sigma_{h}=+1$, asserts $k$ is even. Hence $k=2k^{\prime}$ in (4.27), and transverse crossing (4.12) implies

[TABLE]

for small $a-1>0$. Absence of zero eigenvalues, for $a>1$, together with the strict ordering and pairing of proposition 4.1, (4.6), preserves simplicity of the real eigenvalues and perpetuates (4.29) to all real $a>1$. This proves the proposition. ∎

5 Main result: proof

In this section we return to the original setting

[TABLE]

of our main result, theorem 2.2, with the normalizations (2.9) – (2.11). We also recall the notation $\langle\alpha\rangle,\langle\alpha\rangle_{h},$ of (2.12), (2.17) for the arithmetic and harmonic means of the diagonal elements $-a\alpha_{m}$ of the fast $N$-cycle $\mathbf{A}$, with signs $\sigma,\sigma_{h}=\pm 1$. For the choice

[TABLE]

of the parameter $a$, depending on the feedback sign $\beta=\beta_{1}=\pm 1$, with the remaining off-diagonal elements of $\mathbf{A}$ normalized to $\beta_{2}=\ldots=\beta_{N}=1$, see (2.18), (2.19).

To prove theorem 2.2 we proceed as follows. First we fix the open subset $\mathcal{U}\subseteq J\times\mathbb{R}^{M}$, where we seek global Hopf bifurcation, according to definition 3.2 and corollary 3.4. In lemma 5.1, we then check the crucial assumption (3.11), i.e.

[TABLE]

for the center indices

 

$\Box$ of the Hopf points in $\mathcal{U}$, at $\varepsilon=0.$ An elementary perturbation argument will extend 5.3 to small enough $0<\varepsilon<\varepsilon_{0}$, proving the theorem.

Fix $\varepsilon>0$ small enough. Let $(a,x)\in\mathcal{E}\subseteq J\times\mathbb{R}^{M}$ denote the steady states $f(\varepsilon,a,x)=0$, and distinguish the trivial steady state $f(\varepsilon,a,x^{*})=0$ from the complementary ones:

[TABLE]

Eliminating all nontrivial steady states $\mathcal{E}^{c}$ from further consideration, we define

[TABLE]

as the open background set for global Hopf bifurcation. In other words, the trivial line $\mathcal{E}^{*}$ is the set of steady states in $\mathcal{U}$, and

[TABLE]

is the set of Hopf points in $\mathcal{U}$. Note that $\mathcal{H}$ is finite, by analyticity of the linearization $\mathcal{A}=f_{x}(\varepsilon,a,x^{*})$ in $a$.

5.1 Lemma**.**

Let $\varepsilon=0$. Then the number of Hopf points $(a_{n},x^{*})\in\mathcal{E}^{*}$ for the fast subsystem $\dot{\xi}=\mathbf{A}(a)\xi$ is finite, and

[TABLE]

under any of the assumptions (2.13)–(2.16).

Proof..

By transverse crossings of eigenvalues, in proposition 4.2, we may consider $\underline{a}=0$ in (2.18), and $\bar{a}=1,\infty$ in (2.19), for $\varepsilon=0$, without loss of generality. Let us consider the case $(-1)^{N_{\mathrm{aut}}}=-\beta,\ \bar{a}=\infty$, first, where proposition 4.2 asserts absence of zero eigenvalues of $\mathbf{A}(a)$, for all $a\geq 0$. Then

[TABLE]

all share the same $n$-independent prefactor

[TABLE]

Indeed all strictly unstable dimensions $\mu(a)$ share the same even/odd parity, by absence of any zero eigenvalue. Moreover, the local crossing numbers $\chi(a_{n})$ at $\underline{a}=0<a_{n}<\bar{a}=\infty$ just add up to a net crossing number $\chi\mathrel{\mathop{\mathchar 58\relax}}=\sum\chi(a_{n})$, and therefore

[TABLE]

It remains to show $\mu(\infty)-\mu(0^{+})\neq 0$.

Several subcases arise. Consider the case $\beta=+1$ first. Then our assumption $(-1)^{N_{\mathrm{aut}}}=-\beta$ implies that ${N_{\mathrm{aut}}}$ is odd. Suppose $N\not\equiv 0\ (\mathrm{mod}\ 4)$. Then substitution of (4.21) and (4.18) from proposition 4.3 imply

[TABLE]

On the other hand, the present case $\beta=+1,\ N\not\equiv 0\ (\mathrm{mod}\ 4)$ allows us to invoke the two-line assumption (2.13). Because $N_{\mathrm{aut}}$ is odd, however, the first line of that assumption containing $\sigma_{h}=\pm 1$ cannot hold. The second line of assumption (2.13), which does not contain $\sigma_{h}$, therefore asserts $\mu(\infty)-\mu(0^{+})\neq 0$ in (5.11), as claimed above.

For $\beta=+1,\ N\equiv 0\ (\mathrm{mod}\ 4)$, we similarly substitute (4.22) and (4.18) from proposition 4.3. Invoking the second line of assumption (2.14), without $\sigma_{h}$, then implies

[TABLE]

The cases arising from $\beta=-1$, where ${N_{\mathrm{aut}}}$ is even, are treated analogously, invoking the second lines of assumptions (2.15), (2.16) without $\sigma_{h}$, this time.

It remains to consider the casuistics of $(-1)^{N_{\mathrm{aut}}}=+\beta$. This time, a simple eigenvalue $\lambda_{k}=0$ appears at $a=\bar{a}=1$. Proposition 4.3 guarantees absence of Hopf points, for $a\geq 1$, i.e. $\mu(1+)=\mu(\infty)$. Proposition 4.2 asserts the transverse crossing direction $\mathrm{sign}\,\lambda_{k}^{\prime}(\alpha)=-\sigma_{h}\neq 0$, at $a=1$, i.e. $\mu(1^{+})-\mu(1^{-})=-\sigma_{h}$. Together, this shows

[TABLE]

Due diligence analogous to (5.11) – (5.13), but invoking the first lines of our assumptions (2.13)–(2.16) which contain $\sigma_{h}$ , this time, completes the proof of the lemma. ∎

Proof of theorem 2.2.

Let $(a_{n},x^{*})$ enumerate the finitely many Hopf points of $\mathbf{A}(a)$, at $\varepsilon=0$, ordered such that $0<a_{1}<a_{2}\ldots<a_{\bar{n}}$. Recall that $a=0$ is a Hopf point if, and only if, $N\equiv 1-\beta\ (\mathrm{mod}\ 4)$, by proposition 4.2. Fix any $0<\underline{a}<a_{1}$, in that case, and $\underline{a}=0$, otherwise; see (2.18). Similarly, $\lambda_{k}(a)=0$ occurs, for any $a\geq 0$, if, and only if, $(-1)^{N_{\mathrm{aut}}}=\beta$. By proposition 4.4, we then have $a_{\bar{n}}<1$ and we may fix any $a_{\bar{n}}<\bar{a}<1$, in that case, and $\bar{a}=+\infty$ otherwise; see (2.19). To prove theorem 2.2 we invoke corollary 3.4. In the setting (5.1) – (5.5), it is therefore our only remaining task to show

[TABLE]

for small enough $0<\varepsilon<\varepsilon_{0}$, and for all perturbed Hopf points $(a_{n}(\varepsilon),x^{*})$ of the perturbed matrix family $\mathcal{A}(\varepsilon,a)$ in (5.1). Again, $a_{n}(\varepsilon)\in J=(\underline{a},\bar{a})$ are ordered such that

[TABLE]

At $\varepsilon=0$, and for any $a\in\mathbb{R}$, the matrix $\mathcal{A}$ is block diagonal, with upper left block $\mathbf{A}(a)$, upper right block $\varepsilon\mathbf{B}=0$, lower left block $\varepsilon\mathbf{C}=0$, and lower right block $\varepsilon\mathbf{D}=0$. Standard perturbation theory then asserts $\mathrm{spec}\ \mathcal{A}$ to be given by two disjoint components:

[TABLE]

see for example [Kato80], section II.6. Uniformity of the spectral splitting for $a\rightarrow\infty$ follows from the diagonal limit of $a^{-1}\mathcal{A}(\varepsilon,a)$. Disjointness, for $\varepsilon_{0}$ small enough and uniformly for $a\in J$, follows from proposition 4.1, and the excision of the only zero eigenvalue of $\mathbf{A}$ at $a=1$ in case $(-1)^{N_{\mathrm{aut}}}=\beta$. By our hyperbolicity assumption on $\mathbf{D}$, in theorem 2.2, only the perturbed part $\mathrm{spec}\,\mathcal{A}+o(1)$ contributes any Hopf points to the sum (5.14), and eigenvalues $\lambda_{k}=0$ remain excluded, for $a\in J=(\underline{a},\bar{a})$ and $0<\varepsilon<\varepsilon_{0}$. Note, however, that the specific finite number $\bar{n}(\varepsilon)$ of Hopf points may fluctuate, due to conceivably nontransverse crossings of the Hopf eigenvalues of $\mathbf{A}$ through the imaginary axis, for some $a\in J$ at $\varepsilon=0$. Nevertheless

[TABLE]

allows summation of the crossing numbers $\chi(a_{n}(\varepsilon),x^{*})$, over $n$, to a net crossing number $\chi$, as in (5.10), (5.13), with

[TABLE]

Here the unstable dimensions $\mu$ are evaluated at the fixed boundaries $\underline{a}$ and $\bar{a}$ specified in (2.18), (2.19), where $\mathcal{A}$ is hyperbolic. Since (5.16) implies

[TABLE]

lemma 5.1 establishes claim (5.14), for small enough $0<\varepsilon<\varepsilon_{0}$. This proves our main result, theorem 2.2. ∎

We conclude this section with a few comments on the limitations of our result. Restrictions on the globality trichotomy (3.13) of the connected component $\mathcal{C}$, in corollary 3.4, are caused by our domain $\mathcal{U}=(J\times\mathbb{R}^{M})\setminus\mathcal{E}^{c}$; see (5.5). Indeed the second option of (3.13) calls intersections of $\partial\mathcal{C}$ with

[TABLE]

global. Here we omit $\bar{a}$ in case $\bar{a}=+\infty$, of course. Let $(a_{n},x_{n})$ be a sequence in $\mathcal{C}$ converging to some $(a_{\infty},x_{\infty})\in\partial\mathcal{U}$, with bounded relevant virtual periods $q_{n}\rightarrow q_{\infty}>0$.

Consider the steady state case $(a_{\infty},x_{\infty})\in\mathcal{E}^{c}$, first. Then $(a_{\infty},x_{\infty})$ is another Hopf point, $x_{\infty}\neq x^{*}$, which we had discarded before. Such a Hopf point may occur on another branch of equilibria, like the nontrivial branch bifurcating from the trivial branch $\mathcal{E}^{*}$ at $a=1$, in case $(-1)^{N_{\mathrm{aut}}}=\beta$. Without further assumptions on such nontrivial steady states, this possibility cannot be excluded.

Let us examine the left boundary $a_{\infty}=\underline{a}$ next. In case $N\not\equiv 1-\beta\ (\mathrm{mod}\ 4)$, the left endpoint $\underline{a}=0$ of $\mathcal{E}^{*}$ is not a Hopf point, for any $0<\varepsilon<\varepsilon_{0}$. In fact we could have safely extended our analysis into negative scaling coefficients $a$. The only reason we did not pursue that direction further was our focus on the sign structure of the nonzero diagonal entries $-a\alpha_{m}$ of $\mathbf{A}$; indeed $N_{\mathrm{aut}}$ counts autocatalytic $\alpha_{m}<0$ in our analysis. Reversing all signs of $\alpha_{m}$, and replacing $N_{\mathrm{aut}}$ by $N-N_{\mathrm{aut}}$, the case of negative $a$ becomes a trivial corollary to time reversal, of course, as discussed at the end of section 2.

In case $N\equiv 1-\beta\ (\mathrm{mod}\ 4)$, a Hopf point in $\mathcal{E}^{*}$ of $\mathbf{A}(a)$ occurs at $a=0$, for $\varepsilon=0$. Without further information on the bifurcation direction of the associated Hopf branch of bifurcating periodic solutions $(a,x)$, we cannot make any assertions concerning the sign of $a$, locally, for small $\varepsilon>0$. We therefore eliminated this case by fixing a left boundary $a=\underline{a}>0$ for our domain $\mathcal{U}$, in assumption (2.18).

Similarly, the right boundary $a=\bar{a}<1$ of (2.19), in case $(-1)^{N_{\mathrm{aut}}}=\beta$, eliminated the simple eigenvalue $\lambda_{k}=0$ of $\mathbf{A}$ at $a=1$ from consideration. Indeed suppose the simple eigenvalue $\lambda_{k}(a)$ of order $\langle\alpha\rangle_{h}\cdot(1-a)$, is of the same order as the perturbation $\varepsilon$. Then we may consider the resulting interaction with $\varepsilon\mathbf{D},\ldots$ as a rank-1 perturbation of the $(M-N+1)\times(M-N+1)$ block matrix

[TABLE]

By pole assignment, this may result in arbitrary spectrum of order $\varepsilon$, including multiple steady state bifurcations and Hopf points. Simple planar examples $N=1,\ M=2$ illustrate this. Our choice of $\bar{a}<1$ circumvents such complications.

Finally, our particular choice of the scaling parameter $a$ prevents meaningful results in case $N=2$. Indeed the resulting matrices

[TABLE]

with $\alpha_{2}=(-1)^{N_{\mathrm{aut}}}/\alpha_{1}$, then provide Hopf points $(a,x^{*})$ if, and only if, $N_{\mathrm{aut}}=1,\ |\alpha_{1}|=|\alpha_{2}|=1,\ \beta=-1$, and $|a|<1$. Such an interval of Hopf points violates our condition that Hopf points be isolated. We therefore consider $N\geq 3$, only, and leave the planar case to elementary ODE courses.

6 Four examples

We illustrate theorem 2.2 with four examples, in subsections 6.1–6.4 below: the Oregonator model of the Belousov-Zhabotinsky reaction, Volterra-Lotka population dynamics, the citric acid or Krebs cycle, and a gene regulatory model for mammalian circadian rhythms. Before we address these specific examples, we recall our basic approach in comparison to existing literature, and comment on some advantages, generalizations, and limitations.

Our result is intended as a quick first test to establish the possibility of sustained autonomous oscillations in a given network. Many results are available which exclude oscillations, particularly within the setting (1.5) of mass action kinetics $r_{j}=k_{j}x^{y_{j}}$. We have already mentioned [HJ74, Mie17, Fei19] above. Notably, the results in [Fei19] aim to hold for all positive values of the reaction coefficients $k_{j}$. The results in [HJ74, Mie17], in contrast, usually restrict the reaction rate coefficients. For an example with detailed balance, we recall the Wegscheider relation (1.11). Complex balance is not a remedy. Indeed, the assumed existence of a complex balance equilibrium $x^{*}$ will, in general, impose certain constraints on the reaction rate coefficients: in fact, there are usually more reaction complexes $y_{j},\bar{y}_{j}$ than metabolites $X_{m}$.

We repeat that it is not our concern here, or below, to run anecdotal numerical simulations for one or the other parameter set of reaction coefficients. General results on sustained oscillations usually assert the existence of parameters for Hopf points, typically via a Routh-Hurwitz criterion in general dimension $M$. Even with contemporary methods of computer algebra, and in small dimensions, this remains a formidable task. See for example [GES05, EEetal15] and the references there. Transversality and nonresonance conditions for local Hopf bifurcation are usually left unchecked.

On the surface, we generalize these results in at least two ways. First, our global approach only requires net crossing numbers $\chi$, alias sums of center indices

 

$\Box$ , rather than detailed local analysis. Second, we allow for quite general reaction rate functions $r_{j}=r_{j}(x)$, rather than just mass action kinetics. That much “generality”, however, comes with a twist. A third, and quite substantial, generalization to reversible fast $N$-cycles with positive feedback arises in the framework of Jacobi systems (1.12), on the basis of proposition 4.1 and [FuOl88]. We do not pursue that direction further, here.

Let us address the twist of “generality”, which is directed against mass action. The very setting (2.7) of a fast $N$-cycle $\mathbf{A}$ in the Jacobian $f_{x}=(f_{mm^{\prime}})$ at steady state $x^{*}$ requires the freedom of a decomposition of the partial derivatives $f_{mm^{\prime}}$, alias the partials $r_{jm}=\partial_{x_{m}}r_{j}(x^{*})$, into the fast $N$-cycle $\mathbf{A}$ and the slow remaining partials of order $\varepsilon$, independently from the fixed rates $r_{j}(x^{*})$ themselves which determine the prescribed steady state $x^{*}$. Already Michaelis-Menten kinetics (1.6) provide such freedom of choice:

[TABLE]

Here we may choose $x_{m}^{*}$ as small as we like to guarantee any required range of $r_{jm}$, even for prescribed $r_{j}(x^{*})$. We thus assumed our choice of $\mathbf{A}$, and the slow-fast decomposition (2.7) on the linear level, to be independent from $x^{*}$. See also the sensitivity analysis [BFie18] and the recent substantial generalizations [Vas20], which are based on the same concept. Evidently such independence of $r_{jm}$ from $r_{j}$ fails in the pure mass action case, where all $c_{jm}=0$.

A second caveat concerns our choice of the distinguished bifurcation parameter $a>0$ in our normalization (2.11), $\ a_{m}=a\alpha_{m},\ \prod\alpha_{m}=\pm 1$, of the diagonal entries $a_{m}=-f_{mm}$ of the fast $N$-cycle $\mathbf{A}$. Already for N=2, this scaling prevented a meaningful discussion of $2$-cycles, because the necessary $2\times 2$ Hopf condition $0=\mathrm{tr}\,\mathbf{A}=a(\alpha_{1}+\alpha_{2})$ became invariant under $a$. For general $N$, for example, consider the presence of an invariant stoichiometric subspace:

[TABLE]

for all $j$, and one or several fixed vectors $\mathbf{c}\neq 0$. Indeed, (6.2) implies time invariance of any affine hyperplane $\mathbf{c}^{T}x=\mathrm{const.}$ under the network ODE (1.1). For the $N$-cycle $\mathbf{A}$, this implies $\mathbf{c}^{T}\mathbf{A}=0$, i.e.

[TABLE]

for all $m\,(\mathrm{mod}\ N)$. In particular, the characteristic polynomial (4.8) then reads

[TABLE]

if we assume all $\mathbf{c}_{m}$ are nonzero on the $N$-cycle $\mathbf{A}$. Thus all eigenvalues $\lambda$ simply scale radially outward with $a$, from $\lambda=0$. Such spectral behavior is adverse to Hopf bifurcation. For this formal reason, for example, we do not treat the replicator equation or Eigen’s hypercycle below, which is normalized to the stochastically motivated invariance $x_{1}+\ldots+x_{N}=1$. Of course, we may select other $1$-parameter paths $a_{m}=a_{m}(a)$ in such cases, which are more hospitable towards global Hopf bifurcation as in corollary 3.4. Or else, we may look for fast $N$-cycles, in the present setting, which are only supported on metabolites $m$ for which $\mathbf{c}_{m}=0$, if any.

6.1 Oregonators

The celebrated standard Oregonator [F07] is the simplest, chemically somewhat realistic, model of the Belousov-Zhabotinsky oscillatory reaction mechanism; see [Zha91, Zha07]. In our notation (1.1) the model can be written as

[TABLE]

with mass action rate laws $r_{1},\ldots,r_{5}$ and a stoichiometrically motivated “fudge factor” $c>0$. More generally, we admit arbitrary monotone rate laws $r_{j}$, e.g. of Michaelis-Menten type. For the linearization $\mathcal{A}=f_{x}(x^{*})$ in (2.7) we readily obtain

[TABLE]

Here we use the abbreviation $r^{\prime}_{j}$ for $r_{jm}$, if $r_{j}=r_{j}(x_{m})$ depends on a single metabolite, only. The only feasible $N$-cycle involving $N=3=M$ metabolites is $\mathbf{m}=(3\ 2\ 1)$,

[TABLE]

notably with a strongly autocatalytic step $r_{3}$. Comparison between (6.6) and (2.7) also tells us to consider $\varepsilon\mathrel{\mathop{\mathchar 58\relax}}=r_{21}$ as a small perturbation of the $3$-cycle $\mathbf{A}$ in (2.9), with the normalizations

[TABLE]

Since $N=3\not\equiv 0,2\ (\mathrm{mod}\ 4)$, theorem 2.2 asserts global Hopf bifurcation as follows. If $\beta=+1$, i.e. for $r^{\prime}_{1}>r_{22}$ at steady state $x^{*}$, we are in case (2.13). Therefore $N_{\mathrm{aut}}\in\{0,1\}$ requires $N_{\mathrm{aut}}=1+\sigma_{h}$, which implies $N_{\mathrm{aut}}=0$ and $\sigma_{h}=-1$. These contradictory requirements exclude the case $\beta=+1$ of a positive feedback cycle (6.7), where $r^{\prime}_{1}$ dominates.

For $r^{\prime}_{1}<r_{22}$, i.e. for a negative feedback cycle $\beta=-1$, in contrast, the restrictions (2.15) are satisfied if, and only if $N_{\mathrm{aut}}=0$ or $N_{\mathrm{aut}}=1=-\sigma_{h}$. Specifically this leads to the two cases

[TABLE]

In conclusion, (6.9) implies global Hopf bifurcation for the generalized Oregonator with $a\in(0,\bar{a})$, any small $1-\bar{a}>0$, and for $\varepsilon\mathrel{\mathop{\mathchar 58\relax}}=r_{21}<\varepsilon_{0}(\bar{a})$ small enough.

6.2 Lotka-Volterra networks

In the introduction we have mentioned the planar Lotka system [Lot1920] for oscillating chemical reactions. Going far beyond that classical “predator-prey” system, Volterra [Vol1931] first studied quadratic systems of the general form

[TABLE]

with $m=1,\ldots,M,\ x_{m}>0$, in the context of ecological population dynamics. See [Oli14] for an excellent survey. Usually $a_{mm^{\prime}},a_{m^{\prime}m}>0$ indicate mutually beneficial cooperation or symbiosis between different species $m$ and $m^{\prime}$, whereas $a_{mm^{\prime}},a_{m^{\prime}m}<0$ model mutually toxic competition. The predator-prey case of Lotka is $a_{mm^{\prime}}\cdot a_{m^{\prime}m}<0$. Following standard ecological wisdom, we assume self-inhibition $a_{mm}<0$ to prevent unlimited grow-up or blow-up of single species.

We may rescale any fixed positive equilibrium $x^{*}$ to $x^{*}_{m}=1$, for all $m$, without loss of generality. Then the linearization of (6.10) at $x^{*}$ is given by

[TABLE]

In particular we may examine any (relabeled) feedback $N$-cycle $\mathbf{m}=(1\ldots N)$,

[TABLE]

with normalized diagonal

[TABLE]

i.e. $N_{\mathrm{aut}}=0$. In particular note $\sigma=\mathrm{sign}\,\langle\alpha\rangle=\mathrm{sign}\,\langle 1/\alpha\rangle=\sigma_{h}=+1$.

For positive feedback $N$-cycles $\beta=+1$, conditions (2.13), (2.14) boil down to

[TABLE]

Indeed, suppose $N\not\equiv 0\ (\mathrm{mod}\ 4)$. Then $N_{\mathrm{aut}}=0$ in (2.13) and $\sigma_{h}=+1$ require $2\lfloor N/4\rfloor+1>1$, i.e. $N\geq 5$, for $N\not\equiv 0\ (\mathrm{mod}\ 4)$. For $N\equiv 0\ (\mathrm{mod}\ 4)$, condition (2.14) and $\sigma=\sigma_{h}=+1$ similarly requires $N/4=0$ or $|2N/4-1|>1$, i.e. $N/4\geq 2$. This proves claim (6.14).

For negative feedback $N$-cycles $\beta=-1$, conditions (2.15), (2.16) analogously require

[TABLE]

It remains to specify the parameter region $a\in(\underline{a},\bar{a})$ of global Hopf bifurcation, according to (2.18), (2.19). For case (6.14) we obtain the conditions

[TABLE]

since $N_{\mathrm{aut}}=0$. Similarly, case (6.15) summarizes as

[TABLE]

Assuming other interactions $a_{mm^{\prime}}$ to be of sufficiently small order $\varepsilon$, and hyperbolicity of the diagonal block $\varepsilon\mathbf{D}$ complementary to the $N$-cycle $\mathbf{m}=(1\ldots N)$, theorem 2.2 implies global Hopf bifurcation for parameters $a\in(\underline{a},\bar{a})$ as described in (6.16), (6.17).

6.3 Citric acid cycles

The citric acid cycle (CAC) or Krebs cycle is a central hub of the oxidative energy metabolism in any cell; see for example [BTGS15], chapter 17. Although variants depend on taxonomy, the following 8-cycle of enzymatic Michaelis-Menten reactions is a central feature:

[TABLE]

Here $X_{1}=$ Citrate, $X_{2}=$ Isocitrate, $X_{3}=\alpha$-Ketoglutarate, $X_{4}=$ Succinyl-coenzyme A, $X_{5}=$ Succinate, $X_{6}=$ Fumarate, $X_{7}=$ Malate, and $X_{8}=$ Oxaloacetate. Side reactions and regulatory influences are omitted. Oscillations have been observed, experimentally, in mitochondria extracts of liver and pancreatic cells; see [MacDetal03]. One motivation is to understand oscillations in insulin production.

In absence of self-regulation, the fast monomolecular feedback $8$-cycle (6.18) with rates $r_{m}=r_{m}(x_{m})$ does not provide global Hopf bifurcation. Indeed, linearization of (6.18) at a steady state $x^{*}$ provides the fast cycle $\mathbf{A}$ in (2.9), with $a_{m}=r^{\prime}_{m}=\beta_{m+1}$. This determines the scaling parameter $a$ to be fixed at $a=1$; see (2.11). Moreover $\lambda_{0}=0$ is the eigenvalue with maximal real part, by its positive (left) eigenvector and for positive feedback $\beta=+1$; see propositions 4.1, (4.2)(i). In particular, the steady state $x^{*}$ is linearly stable and Hopf bifurcation is excluded.

Regulatory and self-regulatory controls of the CAC (6.18), however, are biologically essential. Otherwise energy conversion would run high, for no reason and with nowhere to go. In our setting, regulatory feedbacks are the primary focus. Consider an arbitrary $M$-cycle (6.18), with $m\,(\mathrm{mod}\ M)$. Assume, however, that metabolite $X_{N}$ up- or down-regulates reaction $r_{0}\mathrel{\mathop{\mathchar 58\relax}}\ X_{0}\rightarrow X_{1}$, enzymatically, i.e. without any appreciable effect on the mass balance of $X_{N}$ itself. In other words,

[TABLE]

and $r_{m}=r_{m}(x_{m})$ remains monomolecular for all $m\neq 0$. Deviating from the stereotypical monotonicity assumption (1.7) on the partial derivatives $r_{jm}$ we have to admit $r_{0N}<0$ here, to account for the observed inhibitory regulation of reaction $r_{0}$ by metabolite $X_{N}$. Most importantly, we assume the partial derivative $r_{00}$ to be small of order $\varepsilon$, along with all other partial derivatives outside the fast $N$-cycle $\mathbf{m}=(1\ldots N)$ defined by $r^{\prime}_{1},\ldots,r^{\prime}_{N}$, and $r_{0N}$. This provides a fast $N$-cycle matrix $\mathbf{A}$, as in (2.9), given by

[TABLE]

with $m=1,\ldots,N$. Normalization yields

[TABLE]

In particular our previous Lotka-Volterra discussion of section 6.2 applies, with negative feedback $\beta=-1$. We obtain global Hopf bifurcation for any cycle length $N\geq 3$; see (6.15).

Let us return to the CAC case (6.18) of the experiments in [MacDetal03]. The following four regulatory terms are listed there:

[TABLE]

All these regulations act by enzyme inhibition, i.e. $\beta=-1$. The only exception, excitatory self-regulation $\tilde{r}_{3}(x_{3})$ of $x_{3}$ by the input $x_{3}$ itself, can be subsumed into the definition of the rate $r_{3}$ and is therefore omitted. Regulation with $N=1$ has to be considered small, in our setting, because it is anticipatory, on the wrong side of the diagonal of $\mathbf{A}$, along the 8-cycle (6.20). Therefore (6.15) only provides global Hopf bifurcation by the two inhibitory feedbacks $r_{8}(x_{8},x_{4})$ and $r_{5}(x_{5},x_{8})$, separately, on $a\in(0,\infty)$ and under suitable smallness and nondegeneracy conditions for the large number of remaining entries in a full model of the CAC metabolism.

In summary, our result for $N=4$ in (6.22) points at the inhibitory effect of $X_{4}=$ Succinyl-coenzyme A on the Citrate synthase reaction $r_{8}$, which produces $X_{1}=$ Citrate from $X_{8}=$ Oxaloacetate, as a possible regulatory source of the observed oscillations.

Similarly, $N=3$ in (6.22) points at the inhibitory effect of $X_{8}=$ Oxaloacetate on the Succinate dehydrogenase reaction $r_{5}$ from $X_{5}=$ Succinate to $X_{6}=$ Fumarate, as a second possible regulatory cause of oscillations.

6.4 Circadian gene regulation

Gene regulatory mechanisms for circadian rhythms, on the cell level, have received considerable attention over the past decades; see [Lal17] for some references concerning drosophila. A gene regulatory model for cells in the suprachiasmatic nucleus of mammals was developed, among others, by [Miretal09]; see [MFieMKS13] for simulations of periodic orbits in that model. The model involves a total of $M=21$ components, with $8$ gene activities transcripted intro mRNAs, $8$ corresponding proteins, and $5$ heterodimers of proteins. Below we write gene activities in small italics, and proteins in capitals. Except for dimerizations, all reactions are of Michaelis-Menten type (1.6), with numerous enzyme inhibitory feedback cycles. To indicate such inhibition of a reaction $j\mathrel{\mathop{\mathchar 58\relax}}\ X_{m}\rightarrow X_{\bar{m}}$ by a metabolite $X_{m^{\prime}}$ we write $X_{m^{\prime}}\dashv(X_{m}\rightarrow X_{\bar{m}})$. Gene transcription in itself does not lower gene activity; rather we may consider such steps $j\mathrel{\mathop{\mathchar 58\relax}}\ y_{j}\rightarrow\bar{y}_{j}$ as autocatalytic, $\bar{y}_{jm}=y_{jm}\neq 0$, in the language of (1.3), without depleting $x_{m}$. We indicate such steps by arrows $\mapsto$. Finally, all components are subject to linear decay rates.

We do not bother here to write down the complete model network, the ODE model, or any of the more than $150$ rate coefficients, many of them guesswork anyway. Instead we highlight the following cycles:

[TABLE]

[TABLE]

[TABLE]

Here $m,n\in\{1,2\}$ distinguish two variants of the per, cry genes and PER, CRY proteins. Since all components $x_{m}$ are subject to decay, and in absence of strong autocatalysis, we obtain corresponding fast $N$-cycles $\mathbf{A}$ with positive $\alpha_{m}$, in (2.9). In particular $N_{\mathrm{aut}}=0,\ \sigma=\sigma_{h}=1$, with arbitrary $a>0$, and $\beta=\pm 1$ as follows:

[TABLE]

For the positive fast feedback cycle $N=5$ in (6.24), our analysis (6.16) in the Lotka-Volterra section 6.2 asserts global Hopf bifurcation for scaling parameters $a\in(0,\bar{a})$ and any $\bar{a}<1$, under the usual smallness and hyperbolic nondegeneracy conditions. For the negative fast feedback cycles $N=3,5$ in (6.23), (6.25), in contrast, our analysis (6.17) in the same section 6.2 asserts global Hopf bifurcation for the whole interval of scaling parameters $a\in(0,\infty)$.