Reaction network.

A reaction network consists of N species participating in R reactions with given kinetics. In the graph of complexes of a reaction network (also referred to as C-graph), edges represent the reactions r₁, …, r_R and nodes are the sets of species at both sides of the reaction arrows, termed complexes . The number of nodes in the graph, i.e. the number of complexes of the network, is denoted by M.

In the context of cell signaling, basal or constitutive formation of a protein A is encoded as a pseudo reaction from a zero complex ∅ → A, while degradation of the protein is represented by a reaction to the zero complex A → ∅. The translocation of a species from/to another compartment is also represented by a pseudo reaction. For example, if a protein in the cytoplasm A_cyt is translocated into the nucleus the corresponding pseudo-reaction is A_cyt → A_nuc.

Let us consider the example in Fig 1A, where proteins A and B bind to form the complex AB which is further converted into its active form AB*. In addition, proteins A and B are being constitutively expressed (as indicated by pseudo reactions from the complex ∅), and degraded (as indicated by pseudo reactions to the complex ∅).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

A) An example network where (i) protein A and protein B interact to form a complex which is further activated. The open version (ii) includes constitutive (basal) formation and degradation of both proteins A and B through pseudo reactions from/to the zero complex. B) An example for (i) uniterminal and (ii) biterminal networks. C) Signal transduction motif in its (i) closed version and (ii) semi-diffusive version in which the basal formation and degradation of E1, A and E2 is taken into account.

https://doi.org/10.1371/journal.pcbi.1005454.g001

This reaction network contains N = 4 species with concentrations c₁ = [A], c₂ = [B], c₃ = [AB] and c₄ = [AB*] participating in R = 7 reactions, and its associated graph has M = 6 complexes.

Molecularity matrix.

The molecularity matrix of a reaction network is an N × M matrix denoted by Y such that Y_ij is the molecularity of species i in complex j. The rank of the molecularity matrix is denoted by ρ. For the example network and ρ = 4.

Stoichiometric matrix.

The (species) stoichiometric matrix is an N × R matrix denoted by S such that: where β_ij is the molecularity of the species i in the set of products and α_ij is the molecularity of the species i in the set of educts of reaction j. For the example network:

The stoichiometric matrix S can be computed from the molecularity matrix Y by taking into account that every reaction from complex i to complex j has an associated vector Y_j−Y_i where Y_j and Y_i are respectively the j-th and i-th columns of the matrix Y. For example, the reaction vector associated with r₁ is Y₂−Y₁ and the matrix S can be retrieved from the matrix Y as:

The rank of the stoichiometric matrix is denoted by s.

Mass conservation relations.

We consider systems that can exchange energy with the environment (thermodynamically non-isolated systems), and distinguish between open and closed systems based on whether or not they can exchange matter with the environment [41]. Closed networks, often called conservative networks in the CRNT literature [42, 43], do not exchange mass with the exterior. In open reaction networks, the exchange of matter with the environment can appear explicitly (via pseudo-reactions with the zero complex representing the environment) or not (like in reactions of the type A → 2A).

During a reaction mechanism, conservation relations (linear combinations of concentrations that remain constant during the process) might appear. Semi-positive conservation relations (positive linear combinations of species concentrations left invariant by the reactions), represent conserved moieties, i.e. biochemical units that participate without loss of integrity in a reaction mechanism [44–46]. Here we refer to these semi-positive conservation relations as mass conservation relations.

Let B be a N × λ matrix whose columns span the nullspace of S^T, such that S^T B = 0. The matrix spanning the nullspace of S^T is not uniquely determined, and we choose B such that all its entries are non-negative (note that this is always possible provided that each conservation law represents conservation of a chemical unit or moiety). For the networks under study, each conservation relation corresponds to a conserved moiety (networks with other types of conservation relations are not in the scope of this paper). We denote by the conserved moieties of the network. The number of mass conservation relations is given by λ = N − s where N is the number of species and s is the rank of S.

Let us take a reference concentration vector c₀, compute the vector of mass conservation constants σ = B^T c₀, and define: (1)

The so called reaction polyhedron or stoichiometric compatibility class for the reference concentration vector (c₀) consists of all c ≥ 0 such that (2)

A reaction network is closed if there is a strictly positive row vector ϑ with N elements such that ϑS = 0. A network is open if there is no such vector. Note that a network can be open and still have some mass conservation relations.

For the open system in Fig 1A, the rank of the stoichiometric matrix is s = 4. Since the number of species is also N = 4, we have that λ = 0 and, hence, there are no mass conservation laws. However, if we consider only its closed counterpart without the pseudo-reactions, the corresponding stoichiometric matrix reads: we have that λ = 2, and we obtain: implying that there are two mass conservation relations of the form: where A₀ = c₁₀ + c₃₀ + c₄₀ and B₀ = c₂₀ + c₃₀ + c₄₀ (index 0 indicates initial amount) are the (constant) total amounts of species A and B, respectively, in all their forms. Note that the strictly positive vector ϑ = (1, 1, 2, 2) fulfills ϑS = 0, as it corresponds to a closed network.

Throughout, we assume that the dynamics of a reaction network are described by a set of ordinary differential equations (ODEs), which is commonly written in this compact form: (3) as the product of the stoichiometric matrix S and the vector of reaction rates v(c, k). Eq (3) together with an initial condition c_{(t = 0)} = c₀ define an initial value problem.

Multistationarity in biochemical reaction networks.

We say that a biochemical reaction network is multistationary if it has the capacity for multiple positive equilibria, i.e., if there exist k > 0, c₁ > 0 and c₂ > 0 with c₁ ≠ c₂ such that Sv(c₁, k) = Sv(c₂, k) = 0 in Eq (3) and c₁, c₂ belong to the same stoichiometric compatibility class (we use the notation x > 0 to indicate a strictly positive vector, i.e. all entries of x are strictly greater than zero). Note that, without mass conservation, the stoichiometric compatibility class of any positive equilibrium is the positive orthant.

In this work we provide methods to detect multi-steady state behavior arising through saddle-node bifurcations (see S1 Appendix for a detailed explanation on the implications of saddle-nodes and saddle-node bifurcations in multistationarity).

Mass action kinetics.

A reaction network is not described completely until we define the kinetics of the participating reactions. In this work we consider mass action kinetics, for which the rates of the reactions are of the form: (4)

Within this framework the number of parameters of the model, i.e. the number of kinetic rate constants k, is equal to the number of reactions R. Depending on the context, as we will see later, it might be more convenient to denote the reaction rates and kinetic constants using the labels of the source (educt) and product complexes of the reactions. Mass action law is frequently preferred in models of signaling pathways, since Michaelis-Menten and similar kinetics employ assumptions on time- and concentration-scale separation that are rarely fulfilled in signaling [47–49].

Main assumptions.

In what follows, we consider only networks fulfilling the assumptions:

1. A.1. Every reaction is endowed with mass action kinetics. As indicated, this is a common assumption in signaling models, where every process can be ultimately represented by mass action kinetics.
2. A.2. The network admits a strictly positive steady state (where the concentrations of all the species are positive). This assumption is also mild in the context of signaling, due to the inherent reversibility of most protein-protein interactions. In practice, a signaling model can always be modified without loss of generality to avoid zero steady states (assuming very low concentrations instead) such that the assumption is fulfilled.

These assumptions are common to both methods presented in this paper.

Additional assumptions.

This (deficiency-oriented) approach is applicable to any cell signaling network fulfilling A.1 and A.2 under an extra mild assumption concerning the graph of complexes, namely, the so-called linkage classes of the network.

The graph of complexes for a reaction network contains ℓ disconnected subgraphs, also called linkage classes .

Two nodes , are said to be strongly linked if there is a directed path from to and also a directed path from to . A terminal strong linkage class is a maximal set of nodes within a disconnected subgraph such that there is no edge pointing to any other set of nodes that are strongly linked.

We say that a network graph (or for simplicity, a network) is uniterminal if every linkage class in the graph contains only one terminal strong linkage class. For illustrative purposes we depict in Fig 1B a linkage class (i) which contains one terminal strong linkage class (and, therefore, the network is uniterminal), and a linkage class (ii) which contains two terminal strong linkage classes (biterminal). In order to apply the deficiency-oriented approach, we need the following assumption to be fulfilled:

1. A.3. The network graph is uniterminal.

Networks which are non-uniterminal are considered, in general, to be unsuited for the description of real chemical systems [21]. In the context of cell signaling it can be said that the vast majority of the networks are uniterminal.

To illustrate this approach let us consider the signal transduction mechanism presented in Fig 1C (i). The subgraph containing nodes is a closed network in which protein A is transformed into its active form A* by enzyme E1, and deactivated by enzyme E2. Furthermore, protein A activates itself through the formation of an intermediate complex AA*. The concentrations of the species involved are c₁ = [A], c₂ = [E1], c₃ = [E2], c₄ = [A*], c₅ = [AE1], c₆ = [A*E2] and c₇ = [AA*].

Labeling of the complexes and reactions.

We label the complexes with numbers from 1 to M (the labeling is arbitrary) before constructing the matrices of the network. In the example, once we label the complexes as in Fig 1C, the molecularity matrix is: and the stoichiometric matrix reads: (5)

Let us introduce a matrix Λ in which entry (i, j) corresponds to node in linkage class , such that:

The dimensions of the matrix Λ are M × ℓ. In the example (considering only the subgraph corresponding to the closed system) there are M = 9 nodes distributed among ℓ = 3 linkage classes such that:

Within the framework of the deficiency-oriented approach, as indicated in the previous section (Mass action kinetics) it is more convenient to denote the reaction rates according to the educt and product complexes, such that the kinetic constant of the reaction from complex j to complex k is denoted as k_jk. In this way, each complex has an associated mass action monomial of the form: (6)

We can define a vector of monomials ψ(c) = (ψ₁, …, ψ_M)^T. The relationship between the vector containing the reaction rates as defined in Eq (4) and the vector of monomials is given by: where diag(k) is an R-diagonal matrix containing the kinetic constants of the reactions and U ∈ {0, 1}^{R × M} is a matrix that selects the source complexes of the reactions [50].

The vector of mass action monomials ψ(c) for the example network reads:

Model equations.

Within this framework the dynamics Eq (3) of a reaction network can be encoded in an equivalent set of ODEs of the form: (7) where Y is the molecularity matrix, A contains the kinetic rate constants and ψ(c) is the vector of mass action monomials Eq (6). Matrix A is built such that its diagonal elements A_{i, i} contain the negative of the sum of the kinetic rate constants corresponding to the reactions going out of complex , while the off-diagonal elements A_{i, j} contain the kinetic constants of the reactions going from complex to complex . In the example:

Matrix B in Eq (1) can be computed from the stoichiometric matrix or, equivalently for the networks under study, from matrix YA such that (YA)^T B = 0, obtaining:

Therefore, we get the following conservation laws: where A₀, E1₀ and E2₀ are constants which represent, respectively, the total amounts of the species A, E1 and E2, including all their forms. Note that the strictly positive vector ϑ = (1, 1, 1, 1, 2, 2, 2) fulfills ϑS = 0, as it corresponds to a closed network.

Deficiency and equilibrium manifold.

The deficiency of a network is a nonnegative integer number defined as:

According to the Deficiency Zero Theorem [21], if the deficiency of a reaction network is zero and, in addition, the network is weakly reversible (every linkage class is a strong terminal linkage class), then, regardless of the values of the kinetic constants, each stoichiometric compatibility class contains precisely one unique steady state, and this steady state is asymptotically stable. Next we obtain an expression for the locus of equilibria of a reaction network in terms of what we call deficiency parameters.

For the uniterminal networks under study, the deficiency δ coincides with the dimension of the so-called deficiency subspace := Ker(Y) ∩ Im(A). A basis ω₁, …, ω_δ for the deficiency subspace can be computed following [37] and taking into account that: (8)

For the example network a basis is:

Since the vector Aψ is a linear combination of the vectors of the basis of , we have: (9) which contains M − ℓ linearly independent algebraic equations. We refer to the linear coefficients α₁, …, α_δ as deficiency parameters. Substituting the vector ψ by the corresponding mass action monomials, we obtain a set of M − ℓ linearly independent equations in terms of the state vector c, the deficiency parameter vector α and the kinetic parameter vector k, namely: (10) which describes the locus of equilibria of Eq (7). Note that, for fixed k, the dimension of the manifold described by Eq (10) is λ, coinciding with the number of mass conservation laws in the system (since there are N + δ variables and M − ℓ equations, and λ = N − s). The set of equations: (11) where W is given by Eq (1) describes the locus of equilibrium points compatible with the reaction polyhedron fixed by σ. Note that, for fixed k and σ, the system is square with N+δ variables and N+δ equations.

In the example, solving Eq (9) we obtain:

Expressing the mass action monomials through the state variables we get the following equations for the equilibrium manifold Eq (10): (12)

According to the deficiency and the dimension of the equilibrium manifold, we classify the networks into three groups [51]:

1. G.1. Proper networks where N = M − ℓ and, therefore, λ = δ.
2. G.2. Over-dimensioned networks, where N < M − ℓ and, therefore, λ < δ.
3. G.3. Under-dimensioned networks, where N > M − ℓ and, therefore, λ > δ.

For proper networks, once we fix the values of the deficiency parameters α, the equilibrium point is fixed. For over-dimensioned networks, once we fix the values of λ components of the deficiency parameter vector α, the equilibrium point is fixed. Finally, for under-dimensioned networks, we need to fix α and λ − δ components of the state vector such that the equilibrium point is fixed.

The network of the example is under-dimensioned since the dimension of the manifold is λ = 3 and the deficiency is δ = 2. As it can be deduced from Eq (12), as soon as we fix the values of α₁, α₂ and c₂ the equilibrium point is fixed.

Sufficient conditions for a saddle-node in presence of mass conservation.

Once we obtain the equations of the equilibrium manifold and the reaction polyhedron we define the following matrix (in what follows D_x F denotes the Jacobian of F with respect to x): (13) where D_c W = B^T and D_α W = 0. Here it is important to remark that this matrix is always square, of dimensions N + δ, for arbitrary deficiency and manifold dimensions. In the example we obtain a 9 × 9 matrix:

In presence of mass conservation, if there exist c > 0, k > 0 and α such that:

1. C.1.
2. C.2.
3. C.3.

the system Eq (11) has a limit point or saddle-node at the (strictly positive) equilibrium given by c, k. The proof is included in S1 Appendix.

Effective search for multistationarity.

Next we propose a two step algorithm for an efficient and systematic search of multistationarity based on the deficiency criterion. First, we search for the kinetic rate constants and steady state concentrations such that the conditions for a saddle-node (C.1 to C.3) are met. We define the decision vector x for the optimization problem as:

1. x = (k₁, …, k_R, α₁, …, α_λ) in case of proper and over-dimensioned networks
2. x = (k₁, …, k_R, α₁, …, α_δ, c₁, …, c_{λ − δ}) in case of under-dimensioned networks.

The objective function to be minimized (under appropriate constraints) is: such that, at the optimum F_def(x^opt) = 0 the conditions for a saddle-node (C.1 to C.3) are fulfilled. The optimization problem can be formulated as: (14)

The lower and upper bounds for the decision variables are denoted by x_L and x_U, respectively, and the constraint c > 0 is imposed according to assumption A.2. Note that decision variables corresponding to kinetic rate constants have positive lower bounds, while decision variables corresponding to deficiency parameters might take negative values. This optimization problem (14) is non-convex and potentially multi-modal. Consequently, a global optimization algorithm [52] has to be used to search for the optimal solution. Here it is important to remark that, if the network is multistationary, the problem has infinitely many solutions for which F_def(x) = 0.

Finally, in order to verify whether the saddle-node is a saddle-node bifurcation, we start a continuation of equilibrium in forward and backward directions by numerical continuation [53]. In case the saddle node is a saddle node bifurcation, two branches of equilibria emerge (there is multistationarity).

An optimal solution for the example network is found for the set of parameters in Fig 2A. The computation time for the optimization with global scatter search [52] was 40 seconds on an Intel 2.8 GHz Xeon processor. Taken with these parameters, the network is bistable as shown in the bifurcation diagram, where the steady state concentration of the species AA* is plotted against the mass conservation constant σ₃ (also denoted as A₀), which is taken as the bifurcation parameter. The continuation analysis is performed using Cl-Matcont [53].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Bifurcation diagrams corresponding to the example network in Fig 1C in its A) closed and B) open versions.

Tangent bifurcations are indicated by red stars. Stable and unstable equilibria branches are indicated by solid and dashed lines, respectively. Bistability regions are enclosed within the grey rectangles.

https://doi.org/10.1371/journal.pcbi.1005454.g002

Additional assumptions.

This approach is applicable to any cell signaling network fulfilling A.1 and A.2 under an additional assumption ensuring that there is no mass conservation, to be introduced next.

In the original work [23] the absence of mass conservation is guaranteed by considering fully diffusive systems, i.e. systems with input and output flows for every species. In the context of cell signaling this would entail degradation and basal formation of all the species. This assumption does not fit well with signaling pathways, where, for example, single proteins but not signaling complexes are subject to basal formation.

In order to cope with realistic cell signaling models we introduce here what we call semi-diffusive networks. Let us first adopt the notation by [23] and differentiate among true reactions (belonging to the biochemical reaction mechanism), outflow reactions (degradation or flow out of the cell compartment) and feed reactions (basal formation or flow into the cell compartment). In Structural Analysis true reactions are referred to as internal reactions and feed and outflow reactions as exchange reactions [50].

First we calculate the conservation laws of the subsystem consisting only of the true reactions. In semi-diffusive networks, there is degradation of all species, as well as inflow (or basal formation) of one species per conservation law. This species is chosen to be the free form of the corresponding protein participating in the conservation law. If the network is not weakly reversible, we need to ensure that our choice is compatible with A.2. The species in the inflow, i.e., the species which are being constitutively generated, are denoted as key species. Now, we are ready to state the assumption that (together with A.1 and A.2) must be verified:

1. A.4. The reaction network is semi-diffusive.

Importantly, cell signaling pathways usually fit into the category of semi-diffusive networks.

The open network in Fig 1C (ii) is an example of a semi-diffusive network with true reactions r₁, …, r₉, outflow reactions r₁₀, …, r₁₆ corresponding to degradation of each of the species and inflow reactions r₁₇, …, r₁₉ corresponding to basal formation of the key species A, E1 and E2 (one per conservation law).

Labeling of the reactions.

In this approach we enumerate the kinetic rate constants according to the labels of the corresponding reactions. The rate of reaction r_j is denoted by v_j and its corresponding kinetic rate constant by k_j. Thereby, in the example the kinetic constant of the reaction going from complex to complex is denoted in what follows by k₁, since it corresponds to reaction r₁. The reaction rate or flux is denoted by v₁.

Model equations.

The model for the cell signaling pathway can be encoded in a system of ODEs of the form: (15) where K ≥ 0 is the constant inflow term while S_to and v_to(c, k) are the stoichiometric matrix and vector of reaction rates containing only the true and outflow reactions. We build the molecularity matrix Y_to as follows: let Y be the molecularity matrix of the true reactions subsystem (i.e. the subsystem consisting only of the true reactions). First we add to the matrix Y a column corresponding to the zero complex, and then we add the columns corresponding to the species that, being in the outflow, do not appear alone in a complex of the true reactions subsystem. For the example in Fig 1C, the species A, E1, E2 and A* do not appear alone in a complex participating in a true reaction, and the molecularity matrix of the semi-diffusive network reads:

Note that the tenth column (to the right of the vertical line) corresponds to the environment node.

The stoichiometric matrix S_to is built from the columns of the molecularity matrix, incorporating only the true and outflow reactions, i.e. r₁, …, r₁₆. Taking into account the stoichiometric matrix for the true reactions given by Eq (5), we can write the stoichiometric matrix including also the outflow reactions S_to as: where I_{N × N} is the identity matrix containing the source complexes for every species in the network (in this case N = 7).

The inflow vector K contains the rates of basal formation for A, E1 and E2: where v₁₇ = k₁₇, v₁₈ = k₁₈ and v₁₉ = k₁₉ are the (constant) reaction rates of r₁₇, r₁₈ and r₁₉.

Let us build a matrix Y_r with columns corresponding to the source complexes for the true and outflow reactions, namely:

Jacobian and the injectivity property.

When we denote by f(c, k) the right hand side of Eq (15), at equilibrium we have: (16)

The inflow term K is a constant vector, and the Jacobian of the system reads:

Since we assume mass action kinetics, the Jacobian can be easily computed from the network matrices as:

For a strictly positive steady state concentration c > 0 we define: (17) and at steady state we have:

A network is said to be injective if p(c₁, k) = p(c₂, k) ⇒ c₁ = c₂ for all k > 0. If, on the contrary, there exist c₁ > 0, c₂ > 0 and k > 0 with c₁ ≠ c₂ such that p(c₁, k) = p(c₂, k), the network is not injective. It has been demonstrated by [23] that a network is injective if and only if:

Sufficient conditions for a saddle-node in absence of mass conservation.

For a fully-diffusive network, where all species are present in the inflow and outflow, if there exist c > 0, k > 0 such that:

1. C.4.
2. C.5.

the network is multistationary [23]. Importantly, non-injectivity (C.4) alone is not sufficient to ensure multiple positive steady states. The proof can be found in the original paper by Craciun and Feinberg [23].

For a semi-diffusive network, if there are k > 0, c > 0 such that:

1. C.6.
2. C.7.

the system Eq (16) has a saddle-node at the (strictly positive) equilibrium given by c, k. The proof can be found in S1 Appendix.

Effective search for multistationarity.

Next we propose an algorithm for an easy and systematic detection of multistationarity based on the injectivity criterion. First we formulate an optimization problem to detect a saddle-node. For convenience, we consider the fluxes (steady state reaction rates, as defined in Structural Analysis [44]), as decision variables. Let us denote by μ the vector containing the fluxes of the true and outflow reactions and compute:

We consider the following objective function in terms of the vector of fluxes μ:

Redefining expression (17) in terms of the vector of fluxes: the optimization problem can be formulated as: (18) where the vector of fluxes μ is the decision vector and the lower and upper bounds for the decision variables are denoted by μ_L and μ_U, respectively. Note that μ_L > 0 is imposed according to assumption A.2. At the optimum F_inj(μ^opt) = 0, the conditions C.6 and C.7 for a saddle-node are fulfilled.

This optimization problem (18) is potentially multi-modal and consequently a global optimization algorithm [52] is used to search for the solution.

If a vector μ^opt is found such that F_inj(μ^opt) = 0, we can retrieve the exact coordinates of a saddle-node in terms of the concentration vector and kinetic parameters. An equilibrium continuation (in forward and backward directions) is started from this point to analyze the behavior of the equilibrium curve in the vicinity of the saddle-node. If the saddle-node is also a saddle-node bifurcation (turning point), multistationarity appears.

The algorithm finds an optimal vector of fluxes μ^opt for our example network in Fig 1C such that F_inj(μ^opt) = 0 (values provided in the legend) and, therefore, the network has a saddle-node. The computation time for the optimization with global scatter search [52] was 80 seconds on an Intel 2.8 GHz Xeon processor. Here we set the steady state concentrations to 1, such that k_i = μ_i for i = 1, …, 16 and k₁₇ = p₁, k₁₈ = p₂, k₁₉ = p₃. In this way, the point c = 1 is a limit point bifurcation where a real eigenvalue crosses the imaginary axis. Starting from this point we compute the bifurcation diagram with respect to k₁₇, which is the kinetic rate constant corresponding to basal formation of protein A. As it is shown in Fig 2B, there is a region of bistability enclosed by two limit point bifurcations.

Ternary complex formation network.

First, we consider the mechanism of ternary complex formation as depicted in Fig 4B. considering only the reactions indicated by solid arrows, i.e. those with no inflow/degradation of any of the involved species. The network is reversible (and therefore uniterminal) with three mass conservation laws corresponding to interferon (here denoted by I), the receptor subunit IFNAR1 (R1) and the receptor subunit IFNAR2 (R2) in all their forms. Therefore, we can apply the deficiency-oriented approach to find parameter vectors compatible with bistable behavior (the deficiency of the network is δ = 1). No parameter vector was found leading to multistability. This result is compatible with the Deficiency One Algorithm [22], which precludes multistationarity for this network.

Ternary complex formation network with IFN in excess.

When interferon is in excess over the rest of the species, its concentration can be assumed to be constant, and be embedded into the corresponding kinetic rate. Under this assumption, the resultant network, depicted in Fig 4C (only the reactions indicated by solid arrows) is again of deficiency δ = 1. Applying the deficiency-oriented approach we do not find any parameter vector for which the network is bistable. The result is again compatible with the outcome of the Deficiency One Algorithm [22] which precludes multiple steady states for any value of the parameters.

Semi-diffusive network with constant IFN inflow.

Next, we consider the mechanism of ternary complex formation taking into account the degradation of all the involved species, the basal expression of the receptor subunits IFNAR1 and IFNAR2, and a constant inflow of interferon. The resultant network is of deficiency δ = 2. This network has no conservation laws and fulfills the requirements to apply the injectivity-oriented approach. No parameters were found leading to multistationary behavior.

Semi-diffusive network with an excess of IFN.

Finally, basal expression of the receptor subunits IFNAR1 and IFNAR2 was considered together with an excess of interferon, and the degradation of all the species involved. The resultant network is of deficiency δ = 1. Using the same methodology as in the preceding case, no parameters where found leading to multistable behavior.

In summary, none of the topological versions considered for the ternary complex formation was found to show bistable behavior. A test with CoNtRol by Donnell et al. [26] confirms that none of the networks considered in the interferon-receptor complex formation is multistationary. Here it is important to remark that we can not preclude the existence of multiple steady states based on the results obtained by the optimization method. For precluding multistationarity, other injectivity-based methods [26–28, 56] provide conclusive results.

Closed network.

Firstly, we consider the mechanism with no inflow/degradation of any of the species involved. The network is closed (there is a strictly positive row vector in the left kernel of the stoichiometric matrix) and uniterminal. The deficiency is δ = 2. We apply the deficiency-oriented approach and find a saddle-node bifurcation, concluding that the system can exhibit multistationarity. First, we find a zero of the objective function F_def, where we start a continuation of equilibria obtaining two different equilibrium branches.

Semi-diffusive network.

Next, we consider the STAT signalling mechanism taking into account the degradation of all species involved and the basal expression of STAT1 and STAT2, as well as a constant input of active receptor (see Fig 5A). The resultant network is of deficiency δ = 8. The saddle-node conditions based on the injectivity criterion were checked by global optimization methods [52]. Again we find a set of fluxes making the objective function F_inj zero (the system admits a saddle-node). Starting a continuation from the saddle-node we conclude that bistability is preserved in the open system.

All bifurcation diagrams are included in S1 Appendix.