A proof of unlimited multistability for phosphorylation cycles

TL;DR

This paper proves that phosphorylation cycles can have an unlimited number of stable steady states, demonstrating complex multistability behavior in biochemical systems through mathematical analysis.

Contribution

The authors prove a conjecture that phosphorylation cycles can have up to 2[n/2]+1 stable steady states, extending understanding of multistability in biochemical networks.

Findings

Phosphorylation cycles can have arbitrarily many stable steady states.

The paper confirms the maximum number of stable states as 2[n/2]+1.

Mathematical techniques used include geometric singular perturbation theory.

Abstract

The multiple futile cycle is a phosphorylation system in which a molecular substrate might be phosphorylated sequentially n times by means of an enzymatic mechanism. The system has been studied mathematically using reaction network theory and ordinary differential equations. It is known that the system might have at least as many as 2[n/2]+1 steady states (where [x] is the integer part of x) for particular choices of parameters. Furthermore, for the simple and dual futile cycles (n=1,2) the stability of the steady states has been determined in the sense that the only steady state of the simple futile cycle is globally stable, while there exist parameter values for which the dual futile cycle admits two asymptotically stable and one unstable steady state. For general n, evidence that the possible number of asymptotically stable steady states increases with $n$ has been given, which has led to the conjecture that parameter values can be chosen such that [n/2]+1 out of 2[n/2]+1 steady states are asymptotically stable and the remaining steady states are unstable. We prove this conjecture here by first reducing the system to a smaller one, for which we find a choice of parameter values that give rise to a unique steady state with multiplicity 2[n/2]+1. Using arguments from geometric singular perturbation theory, and a detailed analysis of the centre manifold of this steady state, we achieve the desired result.

Equations

$${\mathrm{X}}_{i-1}+\mathrm{E}$$

$$\frac{d{x}_i}{dt}$$

$$\frac{d{y}_{1,i}}{dt}$$

$$\frac{d{y}_{2,i}}{dt}$$

$$\frac{d{x}_{\mathrm{E}}}{dt}$$

$$\frac{d{x}_{\mathrm{F}}}{dt}$$

$$E_{\mathrm{tot}}$$

$$F_{\mathrm{tot}}$$

$$X_{\mathrm{tot}}$$

$$\frac{d{x}_i}{d\tau }$$

$$ϵ\frac{d{y}_{1,i}}{d\tau }$$

$$ϵ\frac{d{y}_{2,i}}{d\tau }$$

$$ϵ\frac{d{x}_{\mathrm{E}}}{d\tau }$$

$$ϵ\frac{d{x}_{\mathrm{F}}}{d\tau }$$

$$X_{\mathrm{tot}}=\sum _{i=0}^n{x}_i+ϵ\sum _{i=1}^n(y_{1,i}+y_{2,i}).$$

$$x^{\prime }$$

$$ϵy^{\prime }$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{ϵ}$$

$$y_{1,i}=K_{1,i}^{-1}{x}_{i-1}{x}_{\mathrm{E}},y_{2,i}=K_{2,i}^{-1}{x}_i{x}_{\mathrm{F}},$$

$$K_{j,i}=\frac{d_{j,i}+k_{j,i}}{a_{j,i}},\text{ for}j=1,2,i=1,\dots ,n.$$

$$E_{\mathrm{tot}}=x_{\mathrm{E}}\left(1+\sum _{\ell =1}^n{K}_{1,\ell }^{-1}{x}_{\ell -1}\right),F_{\mathrm{tot}}=x_{\mathrm{F}}\left(1+\sum _{\ell =1}^n{K}_{2,\ell }^{-1}{x}_{\ell }\right),$$

$$y_{1,i}=\frac{K_{1,i}^{-1}{E}_{\mathrm{tot}}{x}_{i-1}}{1+{\sum }_{\ell =1}^n{K}_{1,\ell }^{-1}{x}_{\ell -1}},y_{2,i}=\frac{K_{2,i}^{-1}{F}_{\mathrm{tot}}{x}_i}{1+{\sum }_{\ell =1}^n{K}_{2,\ell }^{-1}{x}_{\ell }}$$

$$a_{1,i}{x}_{i-1}(E_{\mathrm{tot}}-\sum _{\ell =1}^n{y}_{1,\ell })-(d_{1,i}+k_{1,i})y_{1,i}$$

$$a_{2,i}{x}_i(F_{\mathrm{tot}}-\sum _{\ell =1}^n{y}_{2,\ell })-(d_{2,i}+k_{2,i})y_{2,i}.$$

$$\frac{d{x}_i}{d\tau }$$

$$-\frac{k_{2,i}{K}_{2,i}^{-1}{F}_{\mathrm{tot}}{x}_i}{1+{\sum }_{\ell =1}^n{K}_{2,\ell }^{-1}{x}_{\ell }}+\frac{k_{2,i+1}{K}_{2,i+1}^{-1}{F}_{\mathrm{tot}}{x}_{i+1}}{1+{\sum }_{\ell =1}^n{K}_{2,\ell }^{-1}{x}_{\ell }}$$

$$X_{\mathrm{tot}}=n+1,F_{\mathrm{tot}}=E_{\mathrm{tot}}=1,k_{1,i}=K_{1,i},k_{2,i}=K_{2,i},i=1,\dots ,n.$$

$${\alpha }_i=K_{1,i}^{-1},{\beta }_i=K_{2,i}^{-1},i=1,\dots ,n.$$

$$\frac{d{x}_0}{d\tau }$$

$$\frac{d{x}_i}{d\tau }$$

$$\frac{d{x}_n}{d\tau }$$

$$u=\frac{1+{\sum }_{\ell =1}^n{\beta }_{\ell }{x}_{\ell }}{1+{\sum }_{\ell =1}^n{\alpha }_{\ell }{x}_{\ell -1}}.$$

$$x_i=u^i{x}_0,i=0,\dots ,n.$$

$$n+1=x_0(1+u+\cdots +u^n),\text{hence}{x}_0=\frac{n+1}{1+u+\cdots +u^n}.$$

$$1+\sum _{\ell =1}^n{\beta }_{\ell }{u}^{\ell }{x}_0=u+\sum _{\ell =1}^n{\alpha }_{\ell }{u}^{\ell }{x}_0,$$

$$1+u+\cdots +u^n+\sum _{\ell =1}^n{\beta }_{\ell }{u}^{\ell }(n+1)=u(1+u+\cdots +u^n)+\sum _{\ell =1}^n{\alpha }_{\ell }{u}^{\ell }(n+1).$$

Full text

We prove this conjecture here by first reducing the system to a smaller one, for which we find a choice of parameter values that give rise to a unique steady state with multiplicity $2\lfloor\tfrac{n}{2}\rfloor+1$. Using arguments from geometric singular perturbation theory, and a detailed analysis of the centre manifold of this steady state, we achieve the desired result.

1 Introduction

Post-translational modifications of proteins are ubiquitous in almost all molecular processes at the cellular level. Perhaps the most important post-translational modification process is that of phosphorylation, where no, one or several phosphate groups are attached to specific sites of a protein, thereby creating different protein phosphoforms. It is estimated that more than one third of all proteins in eukaryotes are temporarily phosphorylated [5]. Moreover, the number of phosphorylation sites varies greatly from protein to protein and, for example, exceeds $n=20$ in the case of the tumour suppressor protein p53 [25] and $n=80$ in the case of the tau protein, a microtubule stabiliser [17]. However, the phosphorylation process appears in many different guises. One particularly widespread variant is that of distributive sequential phosphorylation where a protein is phosphorylated one site at a time and in sequential order. Dephosphorylation is the reverse process and removes phosphate groups in the reverse order.

The standard model of distributive sequential phosphorylation consists of $6n$ molecular reactions

[TABLE]

where $n$ denotes the number of phosphorylation sites, ${\rm X}_{i}$, $i=0,\ldots,n$, denotes the phosphoform with $i$ phosphate groups attached, ${\rm E}$ is an enzyme (a kinase) that catalyses the phosphorylation process, ${\rm F}$ is another enzyme (a phosphatase) that catalyses the dephosphorylation process, and ${\rm Y}_{1,i}$ and ${\rm Y}_{2,i}$, $i=1,\ldots,n$, are enzyme-substrate complexes. The model is also known as the multiple futile cycle [27].

Assuming mass-action kinetics [7], the reaction network might be modelled by a system of ordinary differential equations (ODEs) in the concentrations of the substances (the phosphoforms, enzyme-substrate complexes and free enzymes). It leads to a high-dimensional ODE system in $3n+3$ variables (concentrations) and $6n$ constants (one for each reaction). It has been argued, but not proved, that the number of possible stable steady states of the system increases linearly with $n$, thereby allowing the cell a large amount of functional flexibility. The phenomenon has been named *unlimited multistability * [26, 27].

In the case $n=1$, the simple futile cycle, it is known that there is a unique positive steady state for fixed total amounts of substrate (protein) and enzymes, and that this solution is globally asymptotically stable [1]. In the case $n=2$, the dual futile cycle, the conclusion was first reached with the help of simulations that for certain choices of parameters there are three steady states, two of which are stable [21]. It was later proved rigorously that at most three steady states exist and that for certain parameter values there are indeed three [27]. In that paper nothing was proved about the stability of these solutions. Later it was formally proven that for suitable choices of parameters and total amounts of substrate and enzymes these three solutions are hyperbolic, two are asymptotically stable and the third is a saddle [15]. The aim of this paper is to extend these results on stability to general values of $n$.

For general $n$, there can be more than two stable steady states [14]. Evidence that there can be as many as $2\lfloor\frac{n}{2}\rfloor+1$ steady states with $\lfloor\frac{n}{2}\rfloor+1$ of them being stable was presented in [26]. Here $\lfloor x\rfloor$ denotes the floor function, the largest integer smaller than a real number $x$. In more detail, it was shown that steady states are in one to one correspondence with the intersections of two curves in the plane and these curves were plotted numerically. It was also indicated how in a certain formal limit these intersections correspond to the roots of a polynomial in one variable. Starting from these considerations it has been proved analytically that there are parameters for which the number of steady states indicated in [26] exist but their stability was not treated [27].

The number of steady states might be larger than $2\lfloor\frac{n}{2}\rfloor+1$ for certain choices of parameter values. For $n=2$, the bound $2\lfloor\frac{n}{2}\rfloor+1$ is an upper bound, but for $n=3$, there can be as many as five steady states [13].

We prove that there exist parameter values for which the ODE system has $\lfloor\frac{n}{2}\rfloor+1$ stable and $\lfloor\frac{n}{2}\rfloor$ unstable steady states (for fixed total amounts of substrate and enzymes). The line of argument is the following. First we will apply geometric singular perturbation theory (GSPT) to reduce the ODE system in $3n+3$ variables to an ODE system in $n+1$ variables, the phosphoform concentrations. This system is known as the Michaelis-Menten limit (or system) of the original system [15]. For the Michaelis-Menten system we will show that there exists a choice of parameter values for which there is one asymptotically stable steady state with multiplicity $2\lfloor\frac{n}{2}\rfloor+1$. To achieve this, we apply a combination of linear algebra and dynamical systems theory. The main hurdle is to establish asymptotic stability of the steady state as it is not hyperbolic but has a centre manifold of dimension one (for fixed total amount of substrate). The leading non-zero coefficient of the system in the direction of the centre manifold is negative but depends on $n$. Having established this, we again apply perturbation theory to conclude the existence of parameter values for which there are $2\lfloor\frac{n}{2}\rfloor+1$ steady states of which $\lfloor\frac{n}{2}\rfloor+1$ are asymptotically stable and $\lfloor\frac{n}{2}\rfloor$ are unstable. With this in place, we finally lift the steady states to the original system while preserving their stability properties.

2 The basic equations

We consider a protein ${\rm X}$ and let ${\rm X}_{i}$ be the protein with $i=0,\ldots,n$ phosphate groups attached to it. The kinase which catalyses the phosphorylation of ${\rm X}$ is denoted by ${\rm E}$ and the phosphatase which catalyses the dephosphorylation of ${\rm X}$ by ${\rm F}$. The phosphorylation and dephosphorylation of the protein proceeds through the formation of enzyme-substrate complexes. For $i=1,\dots,n$, we let ${\rm Y}_{1,i}$ denote the enzyme-substrate complex formed by ${\rm X}_{i-1}$ and ${\rm E}$, and similarly by ${\rm Y}_{2,i}$ the enzyme-substrate complex formed by ${\rm X}_{i}$ and ${\rm F}$. The reaction mechanism under consideration consists of the reactions

[TABLE]

where the labels of the reactions indicate the reaction rate constants. These are positive real numbers. Let $x_{i}$ be the concentration of ${\rm X}_{i}$ for $i=0,\ldots,n$, $y_{1,i}$ and $y_{2,i}$ the concentrations of ${\rm Y}_{1,i}$ and ${\rm Y}_{2,i}$, respectively, for $i=1,\ldots,n$, and let $x_{\rm E}$ and $x_{\rm F}$ be the concentrations of ${\rm E}$ and ${\rm F}$, respectively. Under the assumption of mass-action kinetics, the evolution equations become

[TABLE]

where $a_{2,0}=d_{2,0}=k_{1,0}=0$, $a_{1,n+1}=d_{1,n+1}=k_{2,n+1}=0$. The quantities

[TABLE]

are conserved. Here $E_{\rm tot}$, $F_{\rm tot}$ and $X_{\rm tot}$ are known as the total amounts of kinase, phosphatase and substrate, respectively. Note that using the equation for the conserved quantities $E_{\rm tot}$ and $F_{\rm tot}$, we might eliminate the variables $x_{\rm E}$ and $x_{\rm F}$ from the right hand side of equations (1)-(5) and thus the evolution equations for these two variables can be discarded. The relation (8) might be used to eliminate a further variable but this will not be done here since it would destroy the symmetry of the system.

We next introduce rescaled variables that generalise the rescaling used in the case of the dual futile cycle [15]. For a parameter $\epsilon>0$, let $x_{\rm E}=\epsilon\widetilde{x}_{\rm E}$, $x_{\rm F}=\epsilon\widetilde{x}_{\rm F}$, $y_{j,i}=\epsilon\widetilde{y}_{j,i}$ and $\tau=\epsilon t$. Denote the derivative with respect to $\tau$ by a prime and drop the tildes for convenience. This leads to the equations

[TABLE]

Correspondingly, we get rescaled conserved quantities $E_{\rm tot}$ and $F_{\rm tot}$ identical to (6) and (7) (after cancelling $\epsilon$), and the conserved quantity

[TABLE]

The expressions in the rescaled variables remain regular as $\epsilon$ tends to zero, in the sense that the functions on the right hand sides of the equations extend in a $\mathcal{C}^{\infty}$-manner to $\epsilon=0$, and in the limit some of the evolution equations reduce to algebraic equations. This is an example of the standard situation considered in GSPT. (For an introduction to GSPT see [19][Chapter 1].) In general we have a system of equations of the form

[TABLE]

Here the prime denotes the derivative with respect to $\tau$ and $\epsilon$ is a parameter. Introducing a new time variable $t=\tau/\epsilon$ and denoting the derivative with respect to $t$ by a dot, leads to the extended system

[TABLE]

The variable $x$ is generally referred to as the slow variable and $y$ as the fast variable.

Assume that the equation $g(x,y,0)=0$ is equivalent to $y=h_{0}(x)$ for a continuously differentiable function $h_{0}$ so that the zero set of $g$ is a manifold. Then for $\epsilon=0$ the two equations for $x$ and $y$ in (14) are equivalent to the system $x^{\prime}=f(x,h_{0}(x),0)$, which we refer to as the limiting system. Of particular importance are the eigenvalues of the derivative of $g$ with respect to $y$. In this context these are known as transverse eigenvalues.

The central conclusion is that when none of the transverse eigenvalues has zero real part, there exists an invariant manifold for the extended system called the slow manifold with the following properties:

• •

Its restriction to $\epsilon=0$ is the zero set of $g$.

• •

The restriction of the extended system to the invariant manifold is a system which, when written in terms of the time variable $\tau$, depends in a regular manner on $\epsilon$ and agrees with the limiting system for $\epsilon=0$.

The meaning of the word ‘regular’ here is not only that the right hand sides of the equations are functions of the parameter which for $\epsilon=0$ are as differentiable as the original system but also that the equations can be solved for the time derivatives.

Consider now the evolution equations (9)-(13) for the multiple futile cycle with the variables $x_{\rm E}$ and $x_{\rm F}$ eliminated using the rescaled conservation equations, the analogues of (6) and (7). Then the variables $x_{i}$ might be taken as the slow variables $x$ in the general set-up and the variables $y_{j,i}$ as the fast variables $y$. Here, the equation $g(x,y,0)=0$ is equivalent to setting the right hand side of (2) and (3) to zero and eliminating $x_{\rm E}$ and $x_{\rm F}$ using the conservation equations. By setting (2) and (3) to zero we obtain

[TABLE]

where

[TABLE]

Using the conserved quantities $E_{\rm tot}$ and $F_{\rm tot}$, we further have by insertion

[TABLE]

which leads to

[TABLE]

for $i=1,\ldots,n$. These expressions define the function $h_{0}$ above.

The transverse eigenvalues are by definition the eigenvalues of the linearisation $L$ of the right hand side of the evolution equations for the concentrations of the substrate-enzyme complexes with respect to those complexes. That is, the linearisation of the equations for $i=1,\dots,n$,

[TABLE]

The matrix $L$ can be seen as the sum of a diagonal matrix $L_{1}$ and another matrix $L_{2}$. The diagonal elements of $L_{1}$ are $-(d_{1,i}+k_{1,i})$ for the rows corresponding to $y_{1,i}$, $i=1,\ldots,n$, and $-(d_{2,i}+k_{2,i})$ for the rows corresponding to $y_{2,i}$, $i=1,\ldots,n$. The matrix $L_{2}$ is block diagonal with one block corresponding to each of the two enzymes. We only need to consider one of these blocks as the other one can be treated in a strictly analogous manner. Consider the block corresponding to the enzyme ${\rm E}$. All entries of the $i$-th row of the block are equal to $-a_{1,i}x_{i-1}$. The transverse eigenvalues have negative real part because all eigenvalues of $-L$ have positive real part, as shown in the lemma below.

Lemma 2.1**.**

If $M=(m_{ij})$ is a matrix with elements of the form $m_{ij}=a_{ij}+b_{i}$ where $a_{ij}=0$ for all $i\neq j$ and all $a_{ii}$ and $b_{i}$ are positive, then all eigenvalues of $M$ have positive real parts.

The proof is given in Subsection 5.1.

We conclude that GSPT applies to this situation and that there exists a slow manifold. The restriction of the extended system to the slow manifold gives a family of dynamical systems depending regularly on the parameter $\epsilon$. The level sets of the conserved quantity defined by $X_{\rm tot}$ are transverse to the slow manifold and so restricting to a constant value of this conserved quantity also results in a regular parameter-dependent dynamical system, which we call the completely reduced system. Note that the dimension of the restriction of the extended system to the slow manifold is greater by one than that of the completely reduced system.

At this point we need some material from the theory of dynamical systems. Background on this can be found in [22]. Let $\dot{x}=f(x)$ be a system of ODEs on $\mathbb{R}^{n}$, where $f$ is $C^{1}$. Let $x_{0}$ be a steady state of this system, that is, a point where $f(x_{0})=0$. Let $A=Df(x_{0})$ denote the derivative (Jacobian) of the right hand side of the equations at $x_{0}$. Then $\mathbb{R}^{n}$ can be written as the direct sum of three linear subspaces $E_{-}$, $E_{c}$ and $E_{+}$, which are spanned by the real and imaginary parts of the generalised eigenvectors of $A$ corresponding to the eigenvalues whose real parts are negative, zero and positive respectively. These subspaces are called the stable, centre and unstable subspaces at $x_{0}$. If there are no eigenvalues with zero real part, so that $E_{c}$ is trivial, then the point $x_{0}$ is said to be hyperbolic. Suppose now that we have a parameter-dependent system $\dot{x}=f(x,\alpha)$ where $f(x_{0},0)=0$ and $f$ is $\mathcal{C}^{1}$ in its dependence on both $x$ and $\alpha$. The point $x_{0}$ is a steady state of the system for $\alpha=0$. If this steady state is hyperbolic, then it follows that for $\alpha$ close to zero there exists a unique steady state close to $x_{0}$ and that it is hyperbolic. The dimensions of the stable and unstable subspaces of this steady state are independent of $\alpha$ for $\alpha$ small.

Returning to our concrete example, if the limiting system has $k$ hyperbolic steady states, then so does the system defined by applying all three conservation laws to the original system for $\epsilon$ small (playing the role of $\alpha$ here). The stability type of these steady states is preserved, in the sense that the dimensions of the stable and unstable manifolds are independent of $\epsilon$. In particular, if for one of the steady states of the limiting system, the stable subspace is the whole of $\mathbb{R}^{n}$, then this is also true for the corresponding steady state with $\epsilon>0$. Moreover, this is known to imply that the steady state is asymptotically stable. In general a steady state of the restriction of the system to the slow manifold is a steady state of the extended system and the dimension of its stable manifold in the extended system is the sum of the dimension of its stable manifold in the extended system restricted to the slow manifold and the number of transverse eigenvalues with negative real part.

Using the fact that the right hand sides of (10) and (11) are zero for $\epsilon=0$ and the expressions in (16), it follows that the limiting system, which in this case will also be referred to as the Michaelis-Menten (or MM) system, consists of the equations

[TABLE]

for $i=0,\ldots,n$, where we adopt the convention that the symbols $K_{2,0}^{-1}$, $K_{1,0}^{-1}$, $K_{1,n+1}^{-1}$ and $K_{2,n+1}^{-1}$ are defined to be zero.

The rest of the paper is devoted to proving the following theorems.

Theorem 2.2**.**

There exists a choice of positive parameters (reaction rate constants and $X_{\rm tot}$) such that the Michaelis-Menten system admits $2\lfloor\frac{n}{2}\rfloor+1$ steady states in the linear invariant subspace defined by the total amount $X_{\rm tot}$. Furthermore, relatively to this invariant subspace, $\lfloor\frac{n}{2}\rfloor+1$ of these steady states are asymptotically stable and hyperbolic, and $\lfloor\frac{n}{2}\rfloor$ are unstable and hyperbolic with a one-dimensional unstable manifold.

As a consequence, GSPT allows us to conclude that the original system with evolution equations (1)-(5) also admits a choice of positive parameters such that there are $\lfloor\frac{n}{2}\rfloor+1$ asymptotically stable steady states and $\lfloor\frac{n}{2}\rfloor$ unstable steady states. Indeed, given that the transverse eigenvalues are all negative this means that for $\epsilon$ small there exist steady states of the original system corresponding to the steady states of the MM system. They are hyperbolic and have corresponding stability properties. The sinks remain sinks and the dimension of the unstable manifolds of the other steady states remains one. This completes the proof of the theorem below.

Theorem 2.3**.**

There exists a choice of positive parameters (reaction rate constants and total amounts) such that the $n$-site phosphorylation system with evolution equations (1)-(5) admits $2\lfloor\frac{n}{2}\rfloor+1$ steady states in the linear invariant subspace defined by the total amounts $E_{\rm tot}$, $F_{\rm tot}$ and $X_{\rm tot}$. Furthermore, relatively to this invariant subspace, $\lfloor\frac{n}{2}\rfloor+1$ of these steady states are asymptotically stable and hyperbolic, and $\lfloor\frac{n}{2}\rfloor$ are unstable and hyperbolic with a one-dimensional unstable manifold.

3 Proof of Theorem 2.2

This section is devoted to selecting ‘nice’ parameters such that the MM system has $2\lfloor\tfrac{n}{2}\rfloor+1$ positive steady states with the same $X_{\rm tot}$, of which $\lfloor\frac{n}{2}\rfloor+1$ are asymptotically stable. To achieve this, we first select parameter values such that $(1,\dots,1)\in\mathbb{R}^{n+1}_{>0}$ is the only positive steady state of the MM system and has multiplicity $2\lfloor\tfrac{n}{2}\rfloor+1$. We then study the stability of this steady state. To this end, we need some centre manifold theory and this will now be reviewed. Background on this can be found in [4] and [20]. Consider once again the system $\dot{x}=f(x)$ and a steady state $x_{0}$. Suppose that $f$ is of class $\mathcal{C}^{k}$ for some finite $k\geq 1$. We restrict consideration to a small neighbourhood of $x_{0}$. There exist manifolds $V_{-}$, $V_{c}$ and $V_{+}$ of class $\mathcal{C}^{k}$ passing through $x_{0}$ that are invariant under the flow generated by the differential equations and are tangent to $E_{-}$, $E_{c}$ and $E_{+}$, respectively, at $x_{0}$. They are referred to as stable, centre and unstable manifolds of the system at $x_{0}$. $V_{-}$ and $V_{+}$ are unique but $V_{c}$ is in general not unique. Fortunately, this lack of uniqueness is usually not a problem in applications as will also be illustrated in our case (to be argued later). It may be noted in passing that the analogues of the statements just made are in general not true if $\mathcal{C}^{k}$ with $k$ finite is replaced everywhere by $\mathcal{C}^{\infty}$ in the sense that there might not exist a $\mathcal{C}^{\infty}$ centre manifold [20, Section 5.1].

We show that the Jacobian of the MM system evaluated at the steady state has $n-1$ eigenvalues with negative real part and two zero eigenvalues. Due to the conservation equation for $X_{\rm tot}$, this implies that in the system defined by restriction to a fixed value of $X_{\rm tot}$, this point has a one-dimensional centre manifold. We proceed to study this centre manifold and conclude that the steady state is asymptotically stable. Finally, we use a perturbation argument to ensure that modified parameters can be found for which there are $2\lfloor\frac{n}{2}\rfloor+1$ positive steady states, of which $\lfloor\frac{n}{2}\rfloor+1$ are asymptotically stable and the remaining ones unstable.

In the following we assume $n\geq 2$, and note that the case $n=1$ has already been solved [1]. The outline of this section is as follows:

• §3.1

We find parameter values such that the limiting system has one steady state. Furthermore, we give the eigenvalue structure of the Jacobian evaluated at this steady state.

• §3.2

The centre manifold is studied. We show that the ODE system of the limiting system is attractive towards the steady state along the direction of the centre manifold. Combined with the results of Subsection 3.1, this establishes the asymptotic stability of the steady state.

• §3.3

The parameter values are perturbed to establish the existence of $\lfloor\frac{n}{2}\rfloor+1$ asymptotically stable steady states of the MM system and these are then lifted to the original system.

3.1 A steady state with multiplicity $2\lfloor\tfrac{n}{2}\rfloor+1$

We let

[TABLE]

To simplify the notation, we define

[TABLE]

Given arbitrary positive values of $K_{j,i}$ one can always find positive values of $d_{j,i},a_{j,i},k_{j,i}$ corresponding to $K_{j,i}$, see (15). We will therefore be concerned with choosing $K_{j,i}$.

With these definitions the MM system becomes:

[TABLE]

Recall that the sum of these equations is zero.

We start by reducing the steady state system to a polynomial equation in one variable. We equate the right-hand side of (17) to zero. Let

[TABLE]

By setting the first equation in (17) to zero, we obtain $x_{1}=ux_{0}$. Using this relation and $\frac{dx_{1}}{d\tau}=0$ we obtain $x_{2}=u\,x_{1}=u^{2}x_{0}$. By iteration we obtain

[TABLE]

The assumption on the conserved quantity, $X_{\rm tot}=n+1$, implies

[TABLE]

Using (18) and (19) we find

[TABLE]

which after substitution of the value for $x_{0}$ in (20) and multiplication by $1+u+\dots+u^{n}$ gives

[TABLE]

Finally, rearranging the terms results in the following polynomial equation in one variable:

[TABLE]

We conclude that the positive steady states of the MM system are in one-to-one correspondence with the positive solutions of the univariate polynomial equation in (21). Given a solution $u$ of this equation, then $x_{0},\dots,x_{n}$ are found from (20) and (19).

For any choice of values of $\alpha_{i},\beta_{i}$ such that

[TABLE]

the polynomial on the left-hand side of (21) is $(u-1)^{n+1}$ for $n$ even and $(u-1)^{n}(u+1)$ for $n$ odd. So the only positive root is $u=1$, which in turn implies that $x_{0}=1$, and the steady state is $p=(1,\dots,1)$. This solution has multiplicity $2\lfloor\tfrac{n}{2}\rfloor+1$. Note that positive $\alpha_{i},\beta_{i}$ can always be found such that (22) is satisfied.

We move on to show that $\alpha_{i},\beta_{i}$ additionally can be chosen such that the Jacobian of the MM system evaluated at $p$ has rank $n-1$ (hence the zero eigenvalue has multiplicity two) and the non-zero eigenvalues have negative real part. One zero eigenvalue corresponds to the conserved quantity $X_{\rm tot}$. If the non-zero eigenvalues have negative real part, then the stability of the steady state might be understood from studying the behaviour of the system along the direction corresponding to the second zero eigenvalue. This is done in Section 3.2.

For this purpose we introduce a combinatorial quantity. Define

[TABLE]

and $i=0,\ldots,n$. Note that $\gamma_{n}(0)=\gamma_{n}(n)=-\tfrac{1}{n+1}$ for all $n$. Let $\gamma_{n}=(\gamma_{n}(0),\ldots,\gamma_{n}(n))$.

Furthermore, we recall a standard binomial identity [23]:

[TABLE]

In particular, we have

[TABLE]

Additionally, we will apply Pascal’s recurrence:

[TABLE]

which is valid for non-negative integers, $i,n\geq 0$.

Proposition 3.1**.**

For any real $\alpha_{1}$, define

[TABLE]

Then (22) is satisfied and further,

[TABLE]

Proof.

By definition

[TABLE]

For $i=1,\ldots,n$ the relation $\alpha_{i}=\alpha_{1}-\gamma_{n}(i-1)-\tfrac{1}{n+1}$ holds since $\alpha_{i}=\beta_{i-1}$ for $i\geq 2$ and $\gamma_{n}(0)=\frac{-1}{n+1}$. This gives

[TABLE]

For $n$ even we have

[TABLE]

where (25) is applied. If $n$ is odd, then using the expression of $\gamma_{n}(i)$ (23), we obtain:

[TABLE]

where (25) is used twice. This concludes the proof of (22).

It remains to show that $1+\sum_{i=1}^{n}\alpha_{i}=1+\sum_{i=1}^{n}\beta_{i}=\alpha_{1}n$. The first equality is a consequence of definition of $\alpha_{i}$ and $\alpha_{1}=\beta_{n}$. For the second equality, we have

[TABLE]

If $(n+1)\sum_{i=1}^{n-1}\gamma_{n}(i)=2$ then we are done. To show this, apply (24) with $P(j)=n-2j$ for $n$ odd and $P(j)=n$ for $n$ even (recalling that $n>1$):

[TABLE]

This concludes the proof. ∎

According to Proposition 3.1 and the discussion after (22), by choosing $\alpha_{1}>0$ large enough such that $\beta_{i}>0$ for all $i=1,\dots,n$, then $p=(1,\dots,1)$ is a steady state of the MM system. In the remaining part of the text, we consider $\alpha_{i},\beta_{i}$ chosen such that this is the case. We proceed to study the Jacobian matrix $J\in\mathbb{R}^{(n+1)\times(n+1)}$ of the MM system (17) evaluated at $p$, with this choice of parameters.

We obtain the following partial derivatives of the MM system (17), evaluated at $p$:

[TABLE]

where $\alpha_{n+1}=\beta_{0}=0$. The matrix $J$ might be reformulated in terms of two other matrices:

[TABLE]

The rows $2,\dots,n$ of $\alpha_{1}n\,J_{2}$ are zero because the $j$-th entry is $-\alpha_{j}+\alpha_{j}+\beta_{j-1}-\beta_{j-1}=0$. Using Proposition 3.1, the matrix $\alpha_{1}nJ$ becomes the following symmetric matrix:

[TABLE]

The matrix $\widetilde{J}$ is independent of the choice of $\alpha_{1}$. The rows $2,\dots,n$ of this matrix are clearly linearly independent.

Proposition 3.2**.**

Consider the matrix $\widetilde{J}$ in (26).

• •

$\widetilde{J}$* has rank $n-1$. An orthogonal basis of the kernel of $\widetilde{J}$ is formed by the vector $(1,\dots,1)$ and the vector $v=(v_{0},\dots,v_{n})\in\mathbb{R}^{n+1}$ with*

[TABLE]

• •

$\widetilde{J}$* has $n-1$ real negative eigenvalues, counted with multiplicity, and two zero eigenvalues.*

Proof.

It is straightforward to see that $(1,\dots,1)$ belongs to the kernel of $\widetilde{J}$, because the column sums are zero. This vector is linearly independent of $v$. To show that $v$ also is in the kernel of $\widetilde{J}$ we note that the scalar product of $v$ with the rows $i=2,\dots,n$ is

[TABLE]

For rows with index $1$ and $n+1$, we have

[TABLE]

Hence $v$ is in the kernel of $\widetilde{J}$. Since the rank of $\widetilde{J}$ is at least $n-1$, $(1,\dots,1)$ and $v$ form a basis of the kernel. A simple computation shows that the sum of the entries of $v$ is zero, that is, $(1,\ldots,1)\cdot v=0$:

[TABLE]

Hence the basis is orthogonal.

To study the eigenvalues of $\widetilde{J}$ we proceed as follows. Since $\widetilde{J}$ is real and symmetric, all eigenvalues are real. We already know that zero is an eigenvalue with multiplicity two. By Descartes’ rule of signs, if all coefficients of the characteristic polynomial of $\widetilde{J}$ are non-negative, then the number of positive roots is zero. Hence the remaining $n-1$ roots of the characteristic polynomial (counted with multiplicity) must be negative.

To show that all coefficients of the characteristic polynomial of $\widetilde{J}$ are non-negative, we show that for $i=1,\dots,n-1$, all non-zero principal minors of size $i$ of $\widetilde{J}$ have sign $(-1)^{i}$. Clearly, this holds for the principal minors of size $1$, since the diagonal entries of $\widetilde{J}$ are all negative (recall $n>1$).

Consider the determinant of a matrix of size $i$ of the form

[TABLE]

for the following cases:

• (1)

$a_{j}=2$ for all $j=1,\ldots,i$.

• (2)

$a_{1}=1-\tfrac{1}{n}$ and $a_{j}=2$ for $j>1$.

• (3)

$a_{i}=1-\tfrac{1}{n}$ and $a_{j}=2$ for $j<i$.

For the three cases, define respective vectors: (1) $x=(1,\dots,1)$, (2) $x=(n,n-1,\dots,n-i+1)$, and (3) $x=(n-i+1,\dots,n-1,n)$. Then the components of $-A_{i}x$ are non-negative. Using [12, Theorem 5.4] on the matrix $-A_{i}$ and subsequently using [12, Theorem 5.1 2°], we conclude that the principal minors of $-A_{i}$ are non-negative, hence in particular $\det(A_{i})$ is either zero or has sign $(-1)^{i}$.

For a subset $H$ of $\{1,\dots,n+1\}$ of size $n+1-j$, $j=2,\ldots,n-1$, we consider the submatrix $\widetilde{J}_{H}$ of size $j$ of $\widetilde{J}$ obtained by removing the columns and rows with indices in $H$. If $H$ contains $1$ or $n+1$, then the submatrix $\widetilde{J}_{H}$ is a block matrix with blocks of the form (27). Hence the determinant of $\widetilde{J}_{H}$ is the product of the determinants of the blocks, and it is either zero or has sign $(-1)^{j}$.

If $H$ contains neither $1$ nor $n+1$, then $\widetilde{J}_{H}$ has the form

[TABLE]

where $z_{1}=1$ if $2\not\in H$ and $z_{1}=0$ otherwise, $z_{2}=1$ if $n\not\in H$ and $z_{2}=0$ otherwise, and $B$ is a block matrix of size $j-2$ with blocks of the form (27) with diagonal entries equal to $-2$.

If $z_{1}=z_{2}=0$ then $\det(\widetilde{J}_{H})$ equals $\det(B)$ times the determinant of

[TABLE]

Since the later is positive, the determinant has the desired sign.

If $z_{1}=1$ and $z_{2}=0$ (and analogously for $z_{1}=0,z_{2}=1$), then we expand the determinant along the last column and obtain

[TABLE]

Note in this case we necessarily have $j\geq 3$. Let $B_{1}$ denote the matrix obtained by removing the first column and the first row of $B$. Then $B_{1}$ is a block matrix of size $j-3$ with blocks of the form (27) and hence its determinant, if non-zero, has sign $(-1)^{j-3}$. We then have

[TABLE]

By assumption $n\geq 2$. Further, $\det(B)$ has sign $(-1)^{j-2}$ if non-zero, and $\det(B_{1})$ has sign $(-1)^{j-3}$ if non-zero. Hence, the determinant has sign $(-1)^{j}$ if non-zero.

Finally, assume $z_{1}=1$ and $z_{2}=1$. Then $j\geq 4$. The matrix $B$ has at least two blocks. If $B$ has more than two blocks, say $\ell$ blocks $B_{1},\dots,B_{\ell}$, the determinant of $\widetilde{J}_{H}$ agrees with the product of the determinants of $B_{i}$, $i=2,\dots,\ell-1$, times the determinant of a matrix similar to $\widetilde{J}_{H}$ but with only two blocks $B_{1}$ and $B_{\ell}$. Since the determinant of $B_{i}$, $i=2,\dots,\ell-1$, has the desired sign, it is sufficient to show that the determinant of $\widetilde{J}_{H}$ has sign $(-1)^{j}$ whenever $B$ has exactly two blocks:

[TABLE]

We will show by induction in $j$ that the determinant is $(-1)^{j}\tfrac{n-j}{n}$, which has the desired sign. For the induction basis we need to consider the two smallest cases $j=4,5$, and it is also convenient to check separately the case $j=6$. For $j=4$, the matrix $B$ is

[TABLE]

which has determinant $\tfrac{n-4}{n}$. Similarly, for $j=5,6$, we check that the determinant is $-\tfrac{n-5}{n}$ and $\tfrac{n-6}{n}$ respectively.

Assume now that the statement holds for all $4\leq j^{\prime}\leq j-1$. We will prove the statement for $j^{\prime}=j$, with $j\geq 7$. Let $i,i+1$ be the indices of the last column of the first block and the first column of the second block of $B$, respectively. Since $j\geq 7$ and $i\geq 2$, at least one of the two inequalities hold $i\geq 4$ or $j-i\geq 4$. We assume that $i\geq 4$, meaning the first block of $B$ has at least size $3$, and the other case follows symmetrically. The statement will follow by applying the Laplace expansion of the determinant of $\widetilde{J}_{H}$ along the rows $i-1,i$. For this, we let $\widetilde{J}_{H,\{i-1,i\},\{j_{1},j_{2}\}}$ be the matrix obtained by removing the $(i-1)$-th and $i$-th rows and the $j_{1}$-th and $j_{2}$-th columns of $\widetilde{J}_{H}$. Then

[TABLE]

The matrix $\widetilde{J}_{H,\{i-1,i\},\{i-2,i-1\}}$ has one zero column, hence the first term is zero. By the induction hypothesis, \det\Big{(}\widetilde{J}_{H,\{i-1,i\},\{i-1,i\}}\Big{)}=(-1)^{j-2}\frac{n-j+2}{n}. The $(i-2)$-th column of $\widetilde{J}_{H,\{i-1,i\},\{i-2,i\}}$ has only one nonzero entry, equal to one, in the $(i-2)$-th entry. Removing the $(i-2)$-th column and row of this matrix, we obtain a matrix of the form $\widetilde{J}_{H}$ of size $j-3$. Here we use that $i\geq 4$. These considerations give:

[TABLE]

as claimed. This concludes the proof. ∎

3.2 Study of the centre manifold

The next step towards the proof of Theorem 2.2 is to complete the study of the stability properties of the steady state $p$ by inspecting its centre manifold. Recall the MM system with our choice of parameters given in (17), and the two orthogonal vectors spanning the kernel of the Jacobian of this system evaluated at the steady state $p$ (see Proposition 3.2):

[TABLE]

Recall also that the dynamics of (17) around $p$ is confined to the linear subspace $n+1=x_{0}+\dots+x_{n}$. Let us refer to the system within this subspace as the restricted system. Since $v$ satisfies $e\cdot v=0$, the centre subspace of the steady state with respect to the restricted system is spanned by $v$. This implies that the centre manifold of $p$ is one dimensional and admits a parametrisation (at least locally around $p$) of the form

[TABLE]

where

[TABLE]

and all $h_{i}$ vanish at least as fast as $s^{2}$ near [math]. Let $h(s)=(h_{0}(s),\dots,h_{n}(s))$. Since $x_{0}(s)=1+ns$, we have $s(x)=\frac{x_{0}-1}{n}$ for $x$ on the centre manifold (at least locally around $p$). Although the centre manifold may not be uniquely determined this will not cause a problem. The general theory says that for two different choices of the centre manifold, if the system, and hence the functions $h_{i}$, are $\mathcal{C}^{k}$, then the functions $h_{i}$ for the two choices of the centre manifold agree up to order $k$. Thus we choose a fixed centre manifold and the arguments which follow are independent of that choice. Note that since the system itself is $\mathcal{C}^{\infty}$ there exists a centre manifold of class $\mathcal{C}^{k}$ for any finite $k$. To justify the calculations in the following we just need to choose $k$ sufficiently large.

We investigate the equations of the centre manifold by selecting $n-1$ linearly independent vectors $\nu_{1},\dots,\nu_{n-1}$ in $\langle v,e\rangle^{\perp}\subseteq\mathbb{R}^{n+1}$. Using (28) we obtain the equations $\nu_{i}\cdot x(s)=\nu_{i}\cdot h(s)$, $i=1,\dots,n-1,$ or equally, in terms of the variable $x$ (on the centre manifold),

[TABLE]

If $f$ denotes the right-hand side of the MM system (17), then evaluation at $x=x(\tau)$ and differentiation with respect to time $\tau$ gives

[TABLE]

or simply,

[TABLE]

(Note that $x(\tau)$ and $x(s)$ mean two different things.)

Consider the vector field $f$ multiplied by the denominators of $f$ and expressed in terms of the parameterisation of the centre manifold. It takes the form

[TABLE]

Let $w(s)=(w_{0}(s),\dots,w_{n}(s))$. Then, using (30) and the definition of $w$, we find

[TABLE]

where the dependence on $s$ is suppressed. This expression implies that the linear combination on the left side vanishes at least at one order higher than $w_{0}$, since $h^{\prime}$ vanishes at least at order one.

Our goal is to show that the first non-zero coefficient in the Taylor expansion of $w_{0}$ is negative and of order $n+1$ if $n$ is even and of order $n$ if $n$ is odd. If this is so, then it must be that locally around $p$ the flow of the vector field $w(s)$ is attracted towards $p$ as a negative coefficient of $w_{0}(s)$ implies the absolute value $|s|$ is decreasing towards zero. Specifically, we prove the following theorem, where $w^{(m)}(s)$ denotes the $m$-th order term in the Taylor expansion of $w$.

Theorem 3.3**.**

The first non-vanishing term in the Taylor expansion of $w_{0}$ is negative. Specifically, the first non-vanishing term is

[TABLE]

with $M=n+1$ if $n$ is even and $M=n$ if $n$ is odd.

In order to prove Theorem 3.3, we investigate the order terms $w_{i}^{(m)}(s)$ of $w_{i}(s)$ through a series of technical lemmas. We start by introducing a new quantity. Let

[TABLE]

and note that the zeroth and first order terms of $\Gamma_{n}$ vanish. The proof of the next lemma is in Subsection 5.2.

Lemma 3.4**.**

Based on (28), the following two relations hold

[TABLE]

Using Lemma 3.4 and (28), we have the following expressions:

[TABLE]

where the argument $s$ of $h_{i}$ and $\Gamma_{n}$ is omitted to ease notation. The details of the proof of (3.2) are not given. A key ingredient in the proof of Theorem 3.3 is the following function of $s$:

[TABLE]

where the second equality follows from the trivial fact that $e\cdot w=0$ and the definition of $v$. We further obtain (see Subsection 5.3 for a proof):

[TABLE]

which implies that terms of $\chi_{n}$ of order zero, one and two vanish. We also note that $w^{(0)}(s)=w^{(1)}(s)=0$, which trivially follows from (3.2).

Lemma 3.5**.**

Consider the following statements, where $m\geq 1$ and $\nu_{1},\dots,\nu_{n-1}$ is a basis of $\langle v,e\rangle^{\perp}\subseteq\mathbb{R}^{n+1}$:

1. (i)

$w_{0}^{(m-1)}(s)=0$. 2. (ii)

$\nu_{j}\cdot w^{(m)}(s)=0$* for all $j=1,\ldots,n-1$.* 3. (iii)

$w^{(m)}(s)=K_{m}vs^{m}$, that is, $w_{i}^{(m)}(s)=K_{m}v_{i}s^{m}$ for all $i=0,\ldots,n$ and some $K_{m}\in\mathbb{R}$. 4. (iv)

$\chi_{n}^{(m)}(s)=-K_{m}\frac{n(n+1)(n+2)}{6}s^{m}$* for some $K_{m}\in\mathbb{R}$.*

Then, for $m\geq 1$, the following implications hold:

[TABLE]

Proof.

(i) $\Rightarrow$ (ii) follows from (31). If (ii) holds, since $e\cdot w=0$, we have $w^{(m)}\in\langle\nu_{1},\dots,\nu_{n-1},e\rangle^{\perp}=\langle v\rangle$, which implies (iii). Finally, assume (iii), then, using (34), we have

[TABLE]

This completes the proof. ∎

Let

[TABLE]

where $i=1,\ldots,n-1$, and consider the matrix

[TABLE]

Lemma 3.6**.**

If $w^{(m)}=0$ and $\Gamma_{n}^{(m-1)}=0$, then

[TABLE]

Proof.

First of all note that $e\cdot w^{(m)}=0$, and also $(n+1,n,\dots,2,1)\cdot w^{(m)}=0$. Indeed, $(n+1,n,\dots,2,1)=-(0,\dots,n)+(n+1)e$ and hence, by (35),

[TABLE]

by hypothesis.

The system $w^{(m)}=0$ can be written in matrix form as

[TABLE]

with c^{(m)}=\big{(}c_{0}^{(m)},\dots,c_{n}^{(m)}\big{)}^{t} and

[TABLE]

This matrix has rank $n-1$, since the rows $2,\dots,n$ are linearly independent. Since $e\widehat{A}=(n+1,n,\dots,2,1)\widehat{A}=0$, also $e\cdot c^{(m)}=(n+1,n,\dots,2,1)\cdot c^{(m)}=0$. Hence, by deleting the second and the third row of $\widehat{A}$, moving the first row to be the last, the first column to be the last, dividing the first $n-2$ resulting equations by $n$, and reorganising the vector $(h_{1}^{(m)},\dots,h_{n}^{(m)})$, system (37) can be rewritten as

[TABLE]

A straightforward computation shows that

[TABLE]

Multiplying both sides of equality (38) by $B$ gives the expression in the statement. ∎

Let

[TABLE]

be the $k$-th descending factorial for $k=0,\ldots,n-1$. By definition, $\ell_{[k]}=0$ if $k>\ell$. If $k=0$, then $\ell_{[0]}=1$ for all $\ell\geq 0$. Furthermore, $0_{[0]}=1$ and $0_{[k]}=0$ for $k>0$. We introduce the following vectors in $\mathbb{R}^{n+1}$:

[TABLE]

In particular, $e_{[1]}=(0,1,2,\dots,n)$ and $e_{[2]}=(0,0,2,6,\dots,i(i-1),\dots,n(n-1))$. The proof of the following lemma is given in Subsection 5.4.

Lemma 3.7**.**

For $k=0,\ldots,n-1$, if $n$ is even, and for $k=0,\ldots,n-2$, if $n$ is odd, we have

[TABLE]

Furthermore,

[TABLE]

The above lemma shows that $\gamma_{n}$ is orthogonal to $e_{[k]}$ for all $k=0,\ldots,n-1$ if $n$ is even, and for all $k=0,\dots,n-2$ if $n$ odd.

We note also the following identities for any non-negative integer $k\geq 0$ and $i\geq k$:

[TABLE]

Lemma 3.8**.**

Let $M>0$ and assume that $\Gamma_{n}^{(m)}=0$, for all $0\leq m\leq M$. Then we have

1. (i)

$w^{(m)}=0$* for $0\leq m\leq M+1$,* 2. (ii)

for $2\leq m\leq M+1$,

[TABLE]

where $c_{mj}\in\mathbb{R}$, $j=1,\ldots,m$, are constants such that $c_{mm}=(-1)^{m}\frac{2^{m}}{m!}$.

Proof.

(i) The proof is by induction on $m$. By construction $w^{(0)}=0$ for $m=0$. Now assume that $w^{(m)}=0$ for an index satisfying $0\leq m\leq M$. We show that the statement holds for $m+1$. Since $w^{(m)}=0$, then from Lemma 3.5, $w^{(m+1)}=K_{m+1}vs^{m+1}$ and $\chi_{n}^{(m+1)}=-K_{m+1}\frac{n(n+1)(n+2)}{6}s^{m+1}$. On the other hand, by assumption, $\Gamma_{n}^{(j)}=0$ for $0\leq j\leq M$, in particular for $0\leq j\leq m$, and combining this with equation (35) yields

[TABLE]

Hence $K_{m+1}=0$, and thus $w^{(m+1)}=0$. This shows (i).

(ii) is proven by induction in $m$. Consider $m=2$. The form of $h_{i}^{(2)}$ follows from Lemma 3.6, using that $c_{0}^{(2)}=-2n(n-1)s^{2}$, $c_{i}^{(2)}=4ns^{2}$, $i\neq 0,n$ and $c_{n}^{(2)}=-2n(n-1)s^{2}$, see (36). Indeed, $h_{i}^{(2)}$ equals $ih_{1}^{(2)}$ plus the $(i-1)$-th component of the matrix product $B(\tfrac{c_{3}^{(2)}}{n},\dots,\tfrac{c_{n}^{(2)}}{n},c_{0}^{(2)})^{t}$, which is

[TABLE]

In vector notation, this is $h^{(2)}=h_{1}^{(2)}e_{[1]}+2s^{2}e_{[2]}$. Since $h_{1}^{(2)}$ is of the form $c_{21}s^{2}$ for some $c_{21}\in\mathbb{R}$ we have

[TABLE]

and the claim is true for $m=2$.

Now assume the claim is true for all $m^{\prime}$ such that $2\leq m^{\prime}\leq m\leq M$ and consider $m^{\prime}=m+1$. We start by describing $h^{(m+1)}$. Since $w^{(m+1)}=0$ by (i) and $\Gamma_{n}^{(m)}=0$, we might apply Lemma 3.6. Since $\Gamma_{n}^{(m^{\prime})}=0$ for $0\leq m^{\prime}\leq m$ by assumption, and $m\geq 2$, it holds that

[TABLE]

Then, by Lemma 3.6, $h_{i}^{(m+1)}$ equals $ih_{1}^{(m+1)}$ plus the $(i-1)$-th component of the matrix product B\big{(}\tfrac{c_{3}^{(m+1)}}{n},\dots,\tfrac{c_{n}^{(m+1)}}{n},c_{0}^{(m+1)}\big{)}^{t}. That is, for $i=2,\dots,n$,

[TABLE]

where

[TABLE]

and

[TABLE]

Let $g^{(q)}$ be defined component-wise as $g_{i}^{(q)}=\sum_{\ell=i+1}^{n}h_{\ell}^{(q)}$. Then, the above result might be given in vector notation as:

[TABLE]

Using $\sum_{\ell=0}^{n}h_{\ell}^{(q)}=0$ (29), equation (39), and the induction hypothesis, we note that

[TABLE]

where $c_{q0}=c_{q(q+1)}=0$. Consequently,

[TABLE]

Note that $h_{j}^{(q)}(s)=z_{qj}s^{q}$ for some $z_{qj}\in\mathbb{R}$. Using this expression and the induction hypothesis on (40), we obtain

[TABLE]

This expression shows that $h^{(m)}$ takes the form stated in the lemma. The coefficient of $s^{m+1}e_{[m+1]}$ is, by the induction hypothesis,

[TABLE]

as required. ∎

Lemma 3.9**.**

$\Gamma_{n}^{(m)}=0$* for all $m=0,\dots,n-1$ if $n$ even and for all $m=0,\dots,n-2$ if $n$ odd.*

Proof.

Let $M=n-1$ if $n$ is even and $M=n-2$ if $n$ is odd. The proof is by induction in $m$. By construction $\Gamma_{n}^{(0)}=0$. Assume $\Gamma_{n}^{(m^{\prime})}=0$ for all $0\leq m^{\prime}\leq m<M$, and consider $m+1$. By Lemma 3.8(ii) and the induction hypothesis, the vector of coefficients of $s^{m+1}$ in $h^{(m+1)}$ lives in the vector space spanned by the vectors $e_{[1]},\dots,e_{[m+1]}$. Now, using (32) and Lemma 3.7, we have

[TABLE]

since $m+1\leq M$. ∎

We are now in a situation where we can prove Theorem 3.3.

Proof of Theorem 3.3.

Let $M=n+1$ if $n$ is even and $M=n$ if $n$ is odd. By Lemma 3.9, $\Gamma_{n}^{(m)}=0$ for all $0\leq m\leq M-2$ and hence by Lemma 3.8(i), $w^{(m)}=0$ (and in particular $w_{0}^{(m)}=0$), for all $1\leq m\leq M-1$.

By Lemma 3.5, we have $w_{0}^{(M)}=K_{M}ns^{M}$ and

[TABLE]

Furthermore, by Lemma 3.8(ii), the vector of coefficients of $s^{M-1}$ in $h^{(M-1)}$ lives in the vector space spanned by the vectors $e_{[1]},\dots,e_{[M-1]}$ and the coefficient of $e_{[M-1]}$ is $c_{M-1,M-1}$. Using (32) and Lemma 3.7, we have

[TABLE]

By (35), we have

[TABLE]

where we use that the first summand vanishes. This implies that

[TABLE]

and hence the first non-zero coefficient of the Taylor expansion of $w_{0}$ is $n$ times

[TABLE]

If $n$ is even, $M-1=n$, $(-1)^{n}=1$ , and we have by Lemma 3.7 and Lemma 3.8(ii) that,

[TABLE]

Similarly, if $n$ is odd, then $M-1=n-1$ and

[TABLE]

The statement now follows from $w_{0}^{(M)}=K_{M}ns^{M}$. This concludes the proof of Theorem 3.3. ∎

These computations are supported by computations in Maple for up to $n=15$.

3.3 Concluding the argument

We have shown that the MM system admits a steady state of multiplicity $2\lfloor\tfrac{n}{2}\rfloor+1$ with $n-1$ negative eigenvalues and such that the first non-zero term of the Taylor expansion of $w_{0}(s)$ is negative. For brevity, as in the previous section, we denote this multiplicity by $M$. Next we will perturb the system so that the multiple steady state splits into $M$ steady states of multiplicity one.

To do this note first that given any polynomial of degree $n+1$ whose leading term is $1$ and whose constant term is $-1$ it is possible to choose the coefficients $\alpha_{i}$ and $\beta_{i}$ in (21) so as to reproduce the given polynomial. Moreover this can be done in such a way that if the coefficients in the polynomial are varied the coefficients $\alpha_{i}$ and $\beta_{i}$ depend smoothly on the coefficients of the polynomial. (In what follows the term ’smooth’ is used to mean $\mathcal{C}^{\infty}$.) Next we introduce a specific family of polynomials depending on a parameter $\mu$ by replacing the function $(u-1)^{2k+1}$ used previously by

[TABLE]

This works for $n$ even. We construct a similar family of polynomials if $n$ is odd. We can then choose coefficients $\alpha_{i}(\mu)$ and $\beta_{i}(\mu)$ depending smoothly on $\mu$ to reproduce this family. After that we can choose parameters of the MM system depending smoothly on $\mu$ so as to give rise to these coefficients $\alpha_{i}$ and $\beta_{i}$. In the end we have a family of MM systems depending smoothly on the parameter $\mu$.

For $\mu=0$, we have the multiple steady state which has been studied before. For convenience it will be referred to as the bifurcation point. As $\mu$ is varied the root of the polynomial at $u=1$, which corresponds to the bifurcation point, splits into simple roots. As is common in bifurcation theory we introduce a suspended system by adjoining the equation $\mu^{\prime}=0$ to the given evolution equations. This increases the dimension of the system by one. The bifurcation point is also a steady state of the suspended system. Its centre manifold as a solution of that system is of dimension two and is foliated by invariant curves of constant $\mu$. The invariant curve with $\mu=0$ is the centre manifold of the bifurcation point with respect to the MM system studied previously. The invariant curve for a non-zero value of $\mu$ will be referred to as the ‘perturbed centre manifold’ although it should be emphasised that it itself is not the centre manifold of anything. It is a general property of a centre manifold of any steady state that all other steady states sufficiently close to the original one lie on that centre manifold. Thus for $\mu$ small all steady states of the parameter-dependent MM system for a fixed value of $\mu$ close to the bifurcation point lie on the perturbed centre manifold.

The non-zero roots $u$ of the polynomial are simple. This suggests that the corresponding steady states of the MM system might be hyperbolic, but this is not obvious. It will now be proved that it is in fact true. The restriction of the system to the two-dimensional centre manifold of the bifurcation point with respect to the suspended system can be thought of as a one-dimensional dynamical system depending on the parameter $\mu$. Let us write it in the form $s^{\prime}(\tau,\mu)=q(s,\mu)$, where $q(s,0)$ is the restriction of the vector field $f$ to the centre manifold parametrised with $s$, and ′ denotes derivative with respect to $\tau$. The function $q$ has one zero for $\mu=0$ and $M$ zeroes for $\mu\neq 0$. The aim is to show that $q^{\prime}$ is non-vanishing at each zero of $q$ for $\mu\neq 0$ and $\mu$ sufficiently small. This will be proved by contradiction.

If the statement is false, then there is a sequence $\mu_{j}$ such that $\mu_{j}\to 0$ as $j\to\infty$, and such that the function $q(s,\mu_{j})$ has at least one degenerate root. In particular $q(s,\mu_{j})$ has at least $M+1$ roots counted with multiplicity. Then $q^{\prime}$ has at least $M$ (distinct) roots, all of which must lie between the smallest and the largest positive roots of $q$. Continuing in this fashion, it can be concluded that $q^{(M)}(s,\mu_{j})$ has at least one root, which is between the smallest and the largest root of $q(s,\mu_{j})$. Hence by continuity it follows that $q^{(M)}(0,0)=0$. This contradicts Theorem 3.3, saying that $q^{(M)}(0,0)\not=0$, and so in reality $q^{\prime}$ does not vanish at any zero of $q$ for $\mu\neq 0$. Note that it does not matter whether we argue about the vector field $w$ or $f$ in Theorem 3.3. Since $w=rf$ with $r(s)$ an arbitrarily often differentiable positive function with $r(0)=1$, then $f(s)\approx-as^{M}$ for some $a>0$ near $s=0$.

It follows from this argument that for any $\mu\neq 0$ the steady states of the restriction of the dynamical system to the perturbed centre manifold are hyperbolic. Since the other eigenvalues of the linearisation at the bifurcation are negative, it follows by continuity that at nearby points with $\mu\neq 0$, all eigenvalues have non-zero real parts. More information can be obtained by considering the sign of the function $q^{\prime}$ at those points where it is not zero. At each steady state the sign changes as no zeros are degenerate. The sign can be determined by continuity since in Theorem 3.3 we have established the sign in the case $\mu=0$. Namely, since $M$ is odd and $q^{(M)}(0,0)<0$, locally around zero, $q(s,0)$ is positive for $s<0$ and negative for $s>0$. It follows that within the perturbed centre manifold sinks and sources alternate and the outermost steady states are sinks. Hence if these points, considered as steady states of the full MM system, are ordered in a suitable way, sinks alternate with saddle points whose unstable manifolds are one-dimensional and the outermost ones are sinks. This completes the proof of Theorem 2.2.

4 Discussion

Perhaps it is of relevance to point out differences and similarities to previous work. The proof of Theorem 2.3 is substantially different from the proof given in [15] for the case $n=2$. Generalising the proof of [15] seems non-trivial and impracticable already for $n=4$, the first case that would have given results establishing the existence of more than the two stable steady states known for $n=2$.

The present proof revolves around several key points. First of all, a reduction of the system using time scale separation that produces a one-dimensional centre manifold is performed. A similar procedure has been applied to investigate the stability of steady states in a dynamical system modelling a different biological situation, the Calvin cycle [6]. Secondly, we obtain the desired number of steady states from a single bifurcation point. Since the centre manifold is one-dimensional, this allows us to conclude the stability of the steady states as we unfold the bifurcation point by parameter perturbation.

It is natural to wonder to what extent these techniques could be developed into a general procedure or statement, as unlimited multistationarity has been established for other families of systems [18, 11]. For the specific reduction of the system to be applicable, we need the transverse eigenvalues to have negative real parts. This is true not only for this case but in general for reduction of systems by (so-called) removal of non-interacting species [10, 9, 24, 8]. So this part of the proof could potentially be generalised. However, it seems non-trivial to devise techniques to guarantee and analyse a one-dimensional centre manifold for which a bifurcation point exists.

As mentioned in the introduction, the steady states whose stability has been investigated here are not the most general steady states of the multiple futile cycle, in the sense that there might exist more than $2\lfloor\tfrac{n}{2}\rfloor+1$ for some parameter choices. It was proven in [27] that $2n-1$ is an upper bound of the number of positive steady states. It was conjectured in [13] that this bound is sharp and the conjecture was proved in the case $n=3$. For $n\geq 3$ the Michaelis-Menten system does not have so many steady states and so a reduction to that system cannot be used to prove the conjecture. Perhaps there exists some other rescaling leading to a different limiting system that admits $2n-1$ steady states, which can be lifted to the multiple futile cycle. If this were true, then perhaps the stability analysis of the present paper could also be extended to that case.

The strategy of establishing the occurrence of dynamical features of a reaction network by using a reduction in the sense of GSPT is not limited to the case of steady states and their stability. It has been used to prove the existence of periodic solutions of the MAPK cascade, a phosphorylation system more complicated than the multiple futile cycle [16]. In the set-up of GSPT explained above, and under the assumption that there are no purely imaginary transverse eigenvalues, it is sometimes possible to show that the existence of a periodic solution of the limiting system implies the existence of a periodic solution of the original system. A sufficient condition for this is that the periodic solution of the limiting system is hyperbolic, that is, no eigenvalue of the Poincaré map has modulus one. This strategy has also been suggested in the context of reaction networks [2]. If in addition that solution is stable (which in this case means that the modulus of each eigenvalue of the Poincaré mapping is less than one), then the solution of the original system is also stable. In [16] the strategy could not be employed since hyperbolicity could not be proved. Instead an alternative strategy was used which might be more widely applicable. The existence of periodic solutions is often proved by showing that there is a Hopf bifurcation. It turns out that a Hopf bifurcation in the limiting system implies the presence of a Hopf bifurcation in the original system and hence the existence of periodic solutions of the original system, without any hyperbolicity condition being necessary.

5 Proofs of auxiliary lemmas

5.1 Proof of Lemma 2.1

Recall the statement of the lemma: If $M=(m_{ij})$ is a square matrix such that $m_{ij}=a_{ij}+b_{i}$ where $a_{ij}=0$ for all $i\neq j$ and all $a_{ii}$ and $b_{i}$ are positive, then all eigenvalues of $M$ have positive real part.

To prove this result, consider, for $n\geq 1$, the following $n\times n$ matrix:

[TABLE]

where $a_{i},b_{i}>0$. We want to show that all eigenvalues of $A_{n}$ have positive real part. Given a matrix $B\in\mathbb{R}^{n\times n}$ and two sets $I,J\subseteq\{1,\dots,n\}$ of cardinality $r\leq n$, we denote by $B_{I,J}\in\mathbb{R}^{r\times r}$ the submatrix of $B$ consisting of the rows indexed by $I$ and the columns indexed by $J$. Then

• •

A matrix $B\in\mathbb{R}^{n\times n}$ is a P-matrix if all principal minors, $\det(B_{I,I})$ for $I\subseteq\{1,\dots,n\}$ are positive.

• •

A matrix $B\in\mathbb{R}^{n\times n}$ is sign-symmetric if for every pair sets $I,J\subseteq\{1,\dots,n\}$ of cardinality $r\leq n$, it holds that

[TABLE]

It follows from [3], that if $B$ is a P-matrix and sign-symmetric, then all eigenvalues of $B$ have positive real part. Therefore, it is enough to show that $A_{n}$ is a $P$-matrix and sign-symmetric.

We first show that $A_{n}$ is a P-matrix for all $n\geq 1$. Because non-maximal principal minors of $A_{n}$ are of the same form as $A_{n}$ but of smaller size, it is sufficient to prove that $\det(A_{n})>0$ for all $n\geq 1$. We first subtract the last column of $A_{n}$ from all other columns and obtain a new matrix $A^{\prime}_{n}$:

[TABLE]

Let $r_{i}$ denote the $i$-th row of $A_{n}^{\prime}$. We replace the last row of $A^{\prime}_{n}$, $r_{n}$, by the linear combination

[TABLE]

and obtain the matrix:

[TABLE]

The matrix $A^{\prime\prime}_{n}$ is upper-triangular matrix and the diagonal entries are positive. Therefore, the determinant of $A^{\prime\prime}_{n}$ is positive, and as a consequence $\det(A_{n})$ is also positive. This shows that $A_{n}$ is a P-matrix.

To see that $A_{n}$ is sign-symmetric for all $n\geq 1$, we will show that

[TABLE]

for all $I,J\subseteq\{1,\dots,n\}$ of cardinality $r$. Clearly, we only need to check the case $I\neq J$.

If $I\cap J$ contains strictly less than $r-1$ elements, then $\det((A_{n})_{I,J})=\det((A_{n})_{J,I})=0$. Indeed, assume there are two indices $i,j\in I$ that are not in $J$. Then $(A_{n})_{I,J}$ has rows $(b_{i}\quad b_{i}\quad\cdots\quad b_{i})$ and $(b_{j}\quad b_{j}\quad\cdots\quad b_{j})$ and hence the determinant is zero. By choosing two indices $i,j\in J$ that are not in $I$ we argue that $\det((A_{n})_{J,I})=0$ as well.

Therefore, we only need to check (41) when $I,J$ differ in only one index. Let $I\cap J=\{\ell_{1},\dots,\ell_{r-1}\}$ with $\ell_{1}<\ldots<\ell_{r-1}$ and $i,j,s,k$ such that

[TABLE]

where $\ell_{k}<i<\ell_{k+1}$ and $\ell_{s}<j<\ell_{s+1}$. The matrix $(A_{n})_{I,J}$ might be constructed by taking a matrix of type $A_{r-1}$, built from the data $a_{\ell_{1}},\dots,a_{\ell_{r-1}}$ and $b_{\ell_{1}},\dots,b_{\ell_{r-1}}$, and adding a row $(b_{i}\ \dots\ b_{i})$ after the $k$-th row and a column with entries $b_{\ell_{1}},\dots,b_{\ell_{k}},b_{i},b_{\ell_{k+1}},\dots,b_{\ell_{r-1}}$ after the $s$-th column. By applying the permutation that sends the $k$-th row to the first row and the $s$-th column to the first column, we conclude that $\det(A_{n})_{I,J}$ agrees with $(-1)^{k+s}$ times the determinant of a matrix of the form

[TABLE]

with $(\beta_{1},\dots,\beta_{r})=(b_{i},b_{\ell_{1}},\dots,b_{\ell_{r-1}})$ and $(\alpha_{2},\dots,\alpha_{r-1})=(a_{\ell_{1}},\dots,a_{\ell_{r-1}})$.

Similarly, the matrix $(A_{n})_{J,I}$ might be constructed by taking the matrix $A_{r-1}$ as above, and adding a row $(b_{j}\ \dots\ b_{j})$ after the $s$-th row and a column with entries $b_{\ell_{1}},\dots,b_{\ell_{s}},b_{j},b_{\ell_{s+1}},\dots,b_{\ell_{r-1}}$ after the $k$-th column. Therefore, $\det(A_{n})_{J,I}$ agrees with $(-1)^{k+s}$ times the determinant of a matrix of the form $B_{r}$ above, now with $\beta_{1}=b_{j}$.

Therefore, it is enough to show that the sign of $\det(B_{r})$ does not depend on the values of $\alpha_{i}>0,\beta_{i}>0$. We proceed similarly to the argument given for the first statement. We subtract the first column of $B_{r}$ to all other columns and obtain a new matrix $B^{\prime}_{r}$ equal to:

[TABLE]

The determinant of $B_{r}^{\prime}$ is $\beta_{1}\alpha_{2}\cdot\dots\cdot\alpha_{r-1}$. Therefore the sign of $\det(B_{r}^{\prime})$, and hence of $\det(B_{r})$, is $+1$ and is independent of $\alpha_{i},\beta_{i}$, as desired. This concludes the proof.

5.2 Proof of Lemma 3.4

Proof.

Recall that $\beta_{i}=\alpha_{1}-\gamma_{n}(i)-\tfrac{1}{n+1}$ (Proposition 3.1), that $\gamma_{n}(n)=-\tfrac{1}{n+1}$ and that $\sum_{\ell=1}^{n}\ell\gamma(\ell)=0$ (Lemma 3.7). By definition of $v_{\ell}$ and Proposition 3.1, we have

[TABLE]

Then by definition of $\Gamma_{n}(s)$ and using $\sum_{\ell=1}^{n}h_{\ell}=0$ (29), we have

[TABLE]

Using this information, Proposition 3.1 and that $h_{0}=0$, we compute the two factors:

[TABLE]

∎

5.3 Proof of Equation (35)

We will show that

[TABLE]

Consider the expression of $w_{i}$ in (3.2). We find

[TABLE]

Using $e\cdot h=e\cdot v=0$ and $h_{0}=0$ yields

[TABLE]

Combined with Lemma 3.4, this gives

[TABLE]

5.4 Proof of Lemma 3.7

Proof.

Let $P_{n}(\ell)=n$ for $n$ even, and $P_{n}(\ell)=n-2\ell$ for $n$ odd. Then $\gamma_{n}(\ell)=\tfrac{(-1)^{\ell+1}P_{n}(\ell)}{n(n+1)}\binom{n}{\ell}.$ This implies for $k\geq 0$,

[TABLE]

Since $P_{n}(j+k)$ has degree [math] in $j$ for $n$ even and degree $1$ for $n$ odd, the identity (24) implies $\sum_{\ell=0}^{n}\ell_{[k]}\gamma_{n}(\ell)=0,$ for $n$ even and $n-k>0$, that is, $k<n$, and for $n$ odd and $n-k>1$, that is, $k<n-1$. This shows the first part of the statement.

For $n$ even and $k=n$, the above reduces to

[TABLE]

as in the statement. For $n$ odd and $k=n-1$, it reduces to

[TABLE]

This concludes the proof of the lemma. ∎

Acknowledgements.

EF and CW have been supported by the Independent Research Fund of Denmark. ADR is grateful for the hospitality of the Department of Mathematical Sciences, University of Copenhagen, while an important part of this work was being done.