## Abstract

Dual phosphorylation of proteins is a principal component of intracellular signalling. Bistability is considered an important property of such systems and its origin is not yet completely understood. Theoretical studies have established parameter values for multistationarity and bistability for many types of proteins. However, up to now no formal criterion linking multistationarity and bistability to the parameter values characterizing dual phosphorylation has been established. Deciding whether an unclassified protein has the capacity for bistability, therefore requires careful numerical studies. Here, we present two general algebraic conditions in the form of inequalities. The first employs the catalytic constants, and if satisfied guarantees multistationarity (and hence the potential for bistability). The second involves the catalytic and Michaelis constants, and if satisfied guarantees uniqueness of steady states (and hence absence of bistability). Our method also allows for the direct computation of the total concentration values such that multistationarity occurs. Applying our results yields insights into the emergence of bistability in the ERK–MEK–MKP system that previously required a delicate numerical effort. Our algebraic conditions present a practical way to determine the capacity for bistability and hence will be a useful tool for examining the origin of bistability in many models containing dual phosphorylation.

## 1. Introduction

Dual phosphorylation of proteins is an important process in intracellular signalling. For example, the proteins of the family of mitogen-activated protein kinases (MAP kinases) are commonly activated by phosphorylation at a threonine and a tyrosine residue [1]. Thereby, phosphorylation is catalysed by a second protein, the MAP kinase. Often the reaction follows a sequential distributive mechanism, where a single phosphorylation occurs with each encounter of kinase and MAP kinase [2], as opposed to a processive mechanism, where two or more phosphorylations occur with each encounter [1,2]. If dephosphorylation also follows a sequential distributive mechanism, dual phosphorylation of a generic protein *A* can be schematically represented as in figure 1.

By employing mass action kinetics [3,4] or Michaelis–Menten kinetics [5], it is straightforward to obtain a mathematical model in the form of an ordinary differential equation (ODE) system that describes the dynamics of a dual phosphorylation network as depicted in figure 1. The steady states of this ODE system determine the long-term behaviour of the network. The number and stability of these steady states depend on the system's parameters: reaction rate constants and total concentrations of kinase, phosphatase and substrate.

In this work, we focus on multistationarity (the existence of at least two positive steady states), which is a prerequisite for bistability (the existence of two stable positive steady states). Bistability is regarded as an important property of dual phosphorylation networks, as it would, for example, enable properties like sharp response thresholds and hysteresis [6]. Consequently, there is a growing interest in multistationarity (bistability) in dual phosphorylation, and, more generally, in network motifs frequently encountered in systems biology [7–10]. Parameter values are known such that multistationarity (bistability) occurs in dual phosphorylation systems. For mass action models, the multistationarity region in parameter space can be described implicitly [3].

In enzyme kinetics research, enzymatic reactions like the ones under study here are often characterized by three constants, the Michaelis constant (*K*_{m}-value), the catalytic constant (*k*_{c}) and the equilibrium constant (*K*_{eq}), collectively referred to as *enzymatic constants* [11,12]. Establishing multistationarity (bistability) using the enzymatic constants is not straightforward, because numerical studies (with total concentrations as bifurcation parameter) are currently the only way. However, the ranges of total concentrations are usually very small and thus values leading to multistationarity can be difficult to find—if at all. Moreover, these studies can be inconclusive, as failure to establish multistationarity (bistability) numerically does not imply non-existence for all parameter values. Hence, algebraic conditions on the enzymatic constants guaranteeing multistationarity (bistability) for some values of the total concentrations are desirable, and our main result establishes such a condition. We further establish (different) conditions that guarantee uniqueness of the steady state and hence non-existence of bistability.

We start by summarizing the main results in §2. We then discuss the biochemical reactions associated with sequential, distributive dual phosphorylation in §3, and the mathematical properties of a mass action model derived from figure 1 in §4. We finally state and discuss the main results, theorems 5.1 and 5.2 and corollary 5.1, and their relation to earlier results in the literature in §5. We end the paper with a sketch of the proofs of theorems 5.1 and 5.2 in appendix A; the details are given in the accompanying electronic supplementary material.

## 2. Emergence of multistationarity

In our results, theorems 5.1 and 5.2, we use catalytic and Michaelis constants as introduced in figure 1. We arrive at the following conclusions.

Catalytic constants determine the capacity for multistationarity and bistability. By theorem 5.1, if the catalytic constants *k*_{c,1}, *…* , *k*_{c,4} satisfy
2.1then multistationarity is guaranteed for some values of the total concentrations. As multistationarity is a prerequisite of bistability, we conclude that the *k*_{c}-values enable the emergence of bistability.

*k*_{c}- and *K*_{m}-values guarantee uniqueness of steady states. By theorem 5.2, if *k*_{c,1}, *…* , *k*_{c,4} and *K*_{m,1}, *…* , *K*_{m,4} satisfy
2.2
2.3then, for every value of the total concentrations, there exists a unique positive steady state.

## 3. Modelling dual phosphorylation

The system depicted in figure 1 consists of four phosphorylation events: the phosphorylation of *A* and *A*_{p} by kinase *E*_{1} and the dephosphorylation of *A*_{pp} and *A*_{p} by phosphatase *E*_{2}. The sequential distributive phosphorylation of, for example, unphosphorylated *A* by *E*_{1} involves three elementary reactions: reversible binding of *A* and *E*_{1} followed by an irreversible release of mono-phosphorylated *A*_{p} and *E*_{1} (*AE*_{1} → *A*_{p} + *E*_{1}).

We will denote by [*s*] the concentration of species *s*. If mass action kinetics is used, the reaction rates of the elementary reactions are assumed to be proportional to the product of the reactant concentrations. For example, the rate of the reaction *A* + *E*_{1} *→ AE*_{1} is given by *k*_{1}[*A*][*E*_{1}] with rate constant *k*_{1}. Similarly, if *k*_{2} and *k*_{3} are the rate constants of *AE*_{1} *→ A* + *E*_{1} and *AE*_{1} *→ A _{p}* +

*E*

_{1}, correspondingly, then the rates are

*k*

_{2}[

*AE*

_{1}] and

*k*

_{3}[

*AE*

_{1}], respectively. We use the term

*mass action constants*to denote rate constants coming from mass action kinetics from hereon.

Every phosphorylation event can be described by *K*_{m}-, *k*_{c}- and *K*_{eq}-values. For example, *K*_{m,1}, *k*_{c,1} and *K*_{eq,1} characterize the phosphorylation of *A* by *E*_{1} discussed above. These *enzymatic constants* are in a one-to-one correspondence with the rate constants [12]. If the rate constants *k*_{1}, *k*_{2}, *k*_{3} are given, then
and if *K*_{m,1}, *K*_{eq,1}, *k*_{c,1} are given, then
As *k*_{1}, *k*_{2} are defined and there is a one-to-one correspondence between *k*_{1}, *k*_{2}, *k*_{3} and *K*_{m,1}, *k*_{c,1}, *K*_{eq,1}.

In a similar way, one has mass action constants for the phosphorylation of *A*_{p} by *E*_{1} (phosphorylation event 2), for the dephosphorylation of *A*_{pp} by *E*_{2} (event 3) and for the dephosphorylation of *A*_{p} by *E*_{2} (event 4). As in case of event 1, events 2–4 can also be described by enzymatic constants. All constants describing dual phosphorylation (mass action and enzymatic) are given in table 1. As for event 1, enzymatic constants can be uniquely expressed in terms of enzymatic constants and vice versa. The corresponding expressions are given in table 2.

We collect the mass action constants in the vector *k* = (*k*_{1}, *…* , *k*_{12})^{T} and the enzymatic constants in the vector
3.1We note that table 2 defines a bijection between *k* and *κ* (with inverse ). In what follows, we will use the bijections and to switch between mass action and enzymatic constants.

Conditions containing mass action constants can thus be translated into conditions containing enzymatic constants by means of . And, vice versa, conditions containing enzymatic constants can be translated into conditions containing mass action constants using .

## 4. Dynamical system

As described in [13], the reactions given in table 1 define a dynamical system of the form
4.1In (4.1), is the vector of species concentrations, Γ is a (9 × 12)-matrix of rank 6 called the stoichiometric matrix and *r*(*k*, *x*) is the vector of reaction rate functions. All reactions are assumed to be governed by mass action kinetics where elements *r _{i}* of the rate function vector

*r*(

*k*,

*x*) are given in the electronic supplementary material, together with the list of species concentrations and the stoichiometric matrix Γ. The complete system (4.1) defined by the reactions in table 1 with mass action constants is spelled out in the electronic supplementary material, equations (4a)–(4i). By means of , we obtain the ODE system with enzymatic constants (given in the electronic supplementary material, equation (26)).

The solutions of the dynamical system (4.1) are confined to affine linear subspaces. To see this, let *W* be a full rank matrix such that *W*^{T} Γ = 0 (*W* is given in the electronic supplementary material, equation (6)). As by (4.1), solutions *ϕ*(*t*, *x*_{0}) of (4.1) with initial condition *ϕ*(0) = *x*_{0} satisfy . Thus, *ϕ*(*t*, *x*_{0}) is contained in the set , which is an affine linear subspace.

As the elements of *x* represent concentrations of chemical species, we restrict to the non-negative orthant and define the set as
4.2

where *c* is a constant vector in with elements: the total concentration of the kinase (*c*_{1}), the phosphatase (*c*_{2}) and the protein (*c*_{3}). Hence, equation (4.2) encodes the fact that for given values of the conserved moieties (kinase, phosphatase and protein), the species concentrations *x* are contained in a hyperplane of the concentration space. In the literature on chemical reaction network theory, the set *ω*_{c} is called stoichiometric class. It is straightforward to show that *ω*_{c} is forward invariant, compact and convex (electronic supplementary material, lemma 3.1), properties used in some of the proofs.

Steady states of (4.1) are solutions of the polynomial system Γ *r*(*k*, *x*) = 0. The set of all positive steady states of (4.1), that is, the set , has been studied in [14,15]. In [14], it is shown that *V* is a toric variety and admits a rational parametrization. Here, we use the computer algebra software Mathematica to solve the polynomials given in (26) of the electronic supplementary material for the unknowns *x*_{3}, *x*_{4}, *x*_{5}, *x*_{6}, *x*_{8} and *x*_{9} in terms of the enzymatic constants *κ* and the ‘free’ unknowns *x*_{1}, *x*_{2}, *x*_{7}. We collect the ‘free’ unknowns in the vector . The resulting rational function that parametrizes all positive steady states of (4.1) in terms of *κ* given in (3.1) and *x*_{f} is depicted in figure 2. For later reference, we summarize this fact in the following proposition.

### Proposition 4.1.

*Consider system (4.1) defined by the reactions in table 1 with the enzymatic constants as parameters (i.e. system (26) of the electronic supplementary material). A positive vector is a steady state of this system, if and only if there exist values of the enzymatic constants and values of such that* .

Multistationarity refers to the existence of at least two positive steady states within the same set *ω*_{c} [13] and bistability to the existence of two ‘stable’ steady states in *ω*_{c}. In the light of the previous discussion, the existence of multistationarity (and hence bistability) requires the existence of a positive vector such that the function *χ* in figure 2 intersects the corresponding set *w*_{c} in at least two points. Multistationarity, therefore, depends on the values of the conserved moieties. Concerning the intersection we introduce the following result from [14,15].

### Proposition 4.2.

*Consider system (4.1) as defined in table 1 and let be given. Then, for any choice of (positive) values of the rate constants, system (4.1) has at least one and at most three positive steady states, that is*, .

### Proof.

By [14, remark 4.1], system (4.1) has at least one positive steady state for any choice of (positive) values of the rate constants and total concentrations. By [15, theorems 1 and 3], system (4.1) has at most three steady states. Furthermore, these results guarantee that there exist values of the rate constants and of the total concentrations such that system (4.1) has exactly three steady states.

## 5. Results

In our study of the dual phosphorylation system shown in figure 1 and table 1, we consider an arbitrary but fixed vector of enzymatic constants *κ**. Hence, is also fixed. In the next theorem, we show that multistationarity for (4.1) exists if inequality (2.1) is satisfied.

### Theorem 5.1.

*Consider system (4.1) as defined by the reactions in table 1. If the catalytic constants k _{c,1}, … , k_{c,4}, satisfy inequality (2.1), then there exist values of the total concentrations of kinase, phosphatase and protein, such that system (4.1) has three positive steady states. That is, there exists and such that*

*i*= 1, 2, 3.

In the next corollary, we show how to find values of the total concentrations such that multistationarity and potentially bistability occur. Notably, these total concentrations can be determined directly, without resorting to bifurcation analysis. The corollary employs the Jacobian of the right-hand side of the ODEs given in (26) of the electronic supplementary material evaluated at the parametrization . And, in particular, the coefficient of its characteristic polynomial. (This coefficient is studied in detail in the electronic supplementary material, §§4.3 and 4.4.)

### Corollary 5.1.

*Consider system (4.1) as defined by the reactions in table 1. Every element* *of the set*
5.1*defines a positive vector* *such that there exist steady states* *with* , *i* = 1, 2, 3.

The uniqueness of a steady state of (4.1) is guaranteed if inequalities (2.2) and (2.3) are satisfied.

### Theorem 5.2.

*Consider system (4.1) as defined by the reactions in table 1. If the catalytic constants k*_{c,1}, …, *k*_{c,4} *and the K*_{m}-*values K*_{m,1}, …, *K*_{m,4} *satisfy inequalities (2.2) and (2.3), then, for every value* *of the total concentrations of kinase, phosphatase and protein, system (4.1) has a unique positive steady state in ω*_{c}. *That is, the system* Γ *r*(*k*, *x*) = 0, *x* ∈ *ω*_{c} *has a unique solution for every* .

A sketch of the proofs is given in appendix A; the details are given in the electronic supplementary material.

### 5.1. Multistationarity and bistability for the mitogen-activated protein kinase ERK

To demonstrate our results, we revisit the MAP kinase ERK model with kinase MEK and phosphatase MKP [5]. Values of the *k*_{c}- and *K*_{m}-constants for the four phosphorylation events depicted in figure 1 are given and bistability is established in [5, figs. 3–5].

Our results can be used to investigate the origin of multistationarity and bistability. The region in the space of total concentrations (*c*-space), where multistationarity is possible, is a complicated three-dimensional object. In [5] only a few sample points (together with a two-dimensional section in [5, fig. 5]) are provided. Using the set , it is straightforward to obtain as many *c* points as desired.

Numerical methods can be used to investigate the set . Testing whether multistationarity is possible for a given vector *c** requires confirming the feasibility of . It is possible to test for multistationarity in regions of *c*-space. If we let and define the box then multistationarity is possible for total concentrations *c* within this box, if and only if the system , is feasible.

Following this line of thought, we analyse the set for the ERK–MEK–MKP system. In figure 3*a*, we use Mathematica to approximate the set and the restriction to *x*_{f}-values yielding total concentrations in biologically plausible domains such as 200 nM ≤ *c*_{1}, *c*_{2} ≤ 400 nM and 400 nM ≤ *c*_{3} ≤ 1000 nM. For demonstration purposes, the system is analysed at and bistability is verified numerically in figure 3*b*. If 231.3286 nM ≤ *c*_{1} ≤ 342.6247 nM, then the system has two stable and one unstable steady states. The bistability region is further explored in figure 3*c* where the solid lines represent points where saddle-node bifurcation occurs. At *c*_{3} = 623.0 nM, the values of *c*_{1}, *c*_{2} inside the region guarantee multistationarity. Different values of *c*_{3} will lead to a different region in *c*_{1}–*c*_{2} space.

### 5.2. Controlling the emergence of bistability

Using theorem 5.1, we conclude that while the catalytic constants enable multistationarity, it is the total concentrations that realize it (and consequently bistability). Given *k*_{c}-values satisfying (2.1), multistationarity emerges if the total concentrations *c* are chosen as specified in corollary 5.1 (provided the model given in figure 1 and table 1 is a good enough approximation). As rate constants are much harder to modify than total concentrations (both *in vivo* and *in vitro*), this points to a mechanism for controlling the onset of multistationarity and hence bistability.

### 5.3. Relation to other work

Finally, we compare our results to those in [6–8,16]. In [6,16], it is demonstrated that bistability is possible if *k*_{c,1} < *k*_{c,2} and *k*_{c,3} < *k*_{c,4}. Clearly, these *k*_{c}-values satisfy (2.1). Our results, however, imply that multistationarity and hence potentially bistability are also possible for *k*_{c,1} < *k*_{c,2} and *k*_{c,3} = *k*_{c,4}.

Inequalities (2.1)–(2.3) are given in the Supporting Material of [7,8]. In the Supporting Material of [8], it is argued that inequalities (2.2) and (2.3) imply injectivity of Γ *r*(*k*, *x*), and hence the absence of multistationarity. Moreover, our result shows the existence of a unique steady state for any positive value of *c*. And in [7], it is argued that (2.1) implies existence of values of *c* such that multistationarity is possible. Moreover, our result shows how to find values of *c* such that exactly three steady states exist in the set *ω*_{c}.

Thus, the exact number of steady states has not been reported in the Supporting Material of [7,8]. This is only possible by combining inequalities (2.1)–(2.3) with results from [14] (concerning the rational parametrization of all positive steady states), [15] (concerning the bounds on the number of positive steady states within classes *ω*_{c}) and degree theory. Hence both the existence of three steady states if inequality (2.1) is satisfied and the existence and uniqueness of a steady state if inequalities (2.2) and (2.3) are satisfied are reported here for the first time. Moreover, corollary 5.1 which provides a direct way to obtain *c*-values such that system (4.1) has exactly three positive steady states in the corresponding set *ω*_{c} has not been published before.

## Funding statement

C.C. was financially supported in part from BMBF grant Virtual Liver (FKZ 0315744) and the research focus dynamical systems of the state Saxony-Anhalt.

## Appendix A. A sketch of proof

Here, we present the main points of the proofs of theorems 5.1 and 5.2. Recall the function *f*(*k*, *x*) = Γ *r*(*k*, *x*) from (4.1) as defined by the reactions in table 1 and the set *ω*_{c} given in (4.2). To determine the number of solutions to *f*(*k*, *x*) = 0, *x* ∈ *w*_{c} (i.e. the cardinality of ), we analyse *f*(*k*, *x*) restricted to sets *ω*_{c}. For this purpose, let and be two matrices of full column rank, whose columns are an orthonormal basis of im(Γ) and its orthogonal complement respectively. Introducing the coordinates , we obtain and define the function
A1As , it follows that
A2

Here, *g _{η}* is the restriction of

*f*(

*k*,

*x*) to some stoichiometric class

*ω*

_{c}.

For a given and hence a given we pick any and compute . Now, if and only if and hence if and only if *ξ* is such that . Thus, we define the set which is forward invariant, convex and compact (electronic supplementary material, lemma 3.2).

Recall that we assume that the vectors *κ** and are given. Then, we treat *η* as a parameter and look for the number of solutions of
A3

In the electronic supplementary material, lemma 3.3, we showed the equivalence of the following statements: (i) (ii) with (iii) with and .

Hence, the number of solutions of (A 3) equals the cardinality of the set . By proposition (4.2), we have A4

We use a degree theory argument to determine the number of solutions of (A 3). Thereby, we study the determinant of the Jacobian *G _{η}* of

*g*evaluated at a steady state of (A 3). The determinant det(

_{η}*G*) corresponds to the constant coefficient, denoted by

_{η}*b*

_{6}(

*κ*,

*x*

_{f}), of the characteristic polynomial of the Jacobian of

*f*(

*k*,

*x*) evaluated at

*χ*(

*κ*,

*x*

_{f}) (electronic supplementary material, lemma 4.1).

The (*topological or Brouwer*) *degree* of *g _{η}* denoted by deg(

*g*,

_{η}*Ω*

*, 0) is a non-negative integer number that contains aggregate information about the number of solutions of (A 3). In the electronic supplementary material, lemma 5.4, we show that the degree of*

_{η}*g*can be defined in terms of and the set

_{η}*X*, A5

_{η}And in the electronic supplementary material, corollary 5.6, we show that A6

Finally, we use the electronic supplementary material, proposition 4.2, to complete the proofs of theorems 5.1 and 5.2. By proposition 4.2, if the *k*_{c,i} satisfy (2.1), then there exists such that . Hence, by (A 4)–(A 6), one has for some and system (A 3) has three solutions. Likewise, also by proposition 4.2, if and *K*_{m,i} satisfy (2.2) and (2.3), then for all . Hence, by (A 5)–(A 6), one has for all and system (A 3) has a unique solution. As the proofs of theorems 5.1 and 5.2 are completed.

- Received February 12, 2014.
- Accepted February 25, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.