Chemical Reaction Network Theory elucidates sources of multistability in interferon signaling

Bistability has important implications in signaling pathways, since it indicates a potential cell decision between alternative outcomes. We present two approaches developed in the framework of the Chemical Reaction Network Theory for easy and efficient search of multiple steady state behavior in signaling networks (both with and without mass conservation), and apply them to search for sources of bistability at different levels of the interferon signaling pathway. Different type I interferon subtypes and/or doses are known to elicit differential bioactivities (ranging from antiviral, antiproliferative to immunomodulatory activities). How different signaling outcomes can be generated through the same receptor and activating the same JAK/STAT pathway is still an open question. Here, we detect bistability at the level of early STAT signaling, showing how two different cell outcomes are achieved under or above a threshold in ligand dose or ligand-receptor affinity. This finding could contribute to explain the differential signaling (antiviral vs apoptotic) depending on interferon dose and subtype (α vs β) observed in type I interferons.

+D=FJAH Some denitions and results from bifurcation analysis 1

.1 Saddle-node and saddle-node bifurcation
We follow [1] to dene the concepts of limit point (or saddle-node ) and limit point bifurcation (also known as saddle-node bifurcation , tangent bifurcation , fold bifurcation or turning point ) in the context of bifurcation analysis, and explain their implications for the existence of multiple steady states in nonlinear ODE systems.
Consider the system of ODEs: where u is the state vector, β is a real parameter, and the function f : R n × R → R n is smooth.
Let us dene we say that the system (1.1) is at equilibrium if: The set of points x = (u T , β) T satisfying Eq. (1.2) constitutes the equilibrium manifold in the (u, β)-space. An equilibrium continuation consists of computing a solution set M ⊂ R n+1 of the smooth system (1.2) starting from a given point x 0 ∈ M .
Let Q be the n × (n + 1) matrix: Next we introduce the concept of regular point .
As a consequence of the Implicit Function Theorem the following result holds: Lemma 1 [1] Near any regular point (u * T , β * ) T , Eq. (1.2) denes a solution curve M that passes through (u * T , β * ) T and is locally unique and smooth.
The proof of Lemma 1 can be found in [1].
The following result introduces some important properties of the tangent vector to the solution curve at a regular point: Lemma 2 [1] A tangent vector v to the solution curve M at a regular point The proof of Lemma 2 can be found in [1].
Next, we introduce the concept of limit point (or saddle-node ).
Condition (2) prevents (u * T , β * ) T from being a hysteresis point [3] and can be replaced by other equivalent non degeneracy conditions, including conditions 2b and 2c introduced next. Let us assume that (u * T , β * ) T is a limit point of (1.2) with respect to β (i.e. condition 1 is fullled): Condition 2b (From [4], Saddle-Node Bifurcation Theorem). Suppose that the kernel of the linear transformation D u f (u * , β * ) : R n → R n is spanned by the nonzero vector w ∈ R n . Let D 2 uu f (u * , β * )(w, w) be the second Fréchet derivative of f evaluated at (u * T , β * ) T in the directions given by w and w.
If D 2 uu f (u * , β * )(w, w) ∈ R n is non zero and not in the range of D u f (u * , β * ), the limit point (u * T , β * ) T is a limit point bifurcation.
Condition 2c Let v be the tangent vector to the equilibrium curve at (u T , β) T , if during an equilibrium continuation, the sign of v n+1 (there exists a parametrization of M that makes v n+1 zero) changes as we pass through (u * T , β * ) T , the limit point (u * T , β * ) T is also a limit point bifurcation [5].

Saddle-node bifurcation, multistationarity and multistability
If the point (u * T , β * ) T is a limit point bifurcation for (1.2) with respect to the parameter β the dynamic system (1.1) has multiple steady states , see Fréchet derivative generalizes the notion of gradient to multivariate matrix functions.
If the point (u * T , β * ) T with u * being strictly positive is a limit point bifurcation for (1.2) with respect to the parameter β, the dynamic system (1.1) has multiple positive steady states , i.e. the system is multistationary .
Here it is important to note that situations in which a limit point is not a limit point bifurcation (turning point) are quite exceptional in practice.
Multiple steady states (although not frequently) might not imply multistability (multiple stable steady states).
Starting a continuation of equilibrium from a limit point (or saddle-node) in forward and backward directions [5] we elucidate whether the limit point is a saddle-node bifurcation or a hysteresis point. If the point is a saddle-node bifurcation, two branches of equilibria emerge, one stable and one unstable. In order to check whether a multistationary system is multistable, we continue the equilibrium curve until we nd a second saddle-node bifurcation.

+D=FJAH
Sucient condition for a saddle-node in reaction networks with mass conservation For a reaction network with N species in M complexes, participating in R reactions endowed with mass action kinetics, the dynamics are given by a system of N ODEs of the form: where c ∈ R N is the vector of concentrations of the species involved, A is a M × M matrix containing kinetic constants constructed as indicated in the main text and ψ(c) ∈ R M is the vector of mass action monomials. Matrices Y and S are the molecularity and stoichiometric matrices respectively, also introduced in the main text.
Let us denote by ℓ the number of linkage classes of the graph of complexes, and by D the deciency subspace (of dimension δ) of the network.
In this section we consider reaction networks in presence of mass conservation fullling assumptions in the main text, i.e.: A.1 Mass action kinetics.
A.2 Capacity for a positive equilibria.
The locus of equilibria of (2.1) can be expressed by a set of M − ℓ algebraic equations of the form: where α ∈ R δ and k ∈ R R are the deciency and kinetic parameter vectors, respectively. See the main text for more details on how H(c, α, k) is computed. Detailed computations are included for all the examples in following chapters of this S1 Appendix (4.1, 4.2, 5.1 and 7).
Let s be the rank of the stoichiometric matrix. For networks with conservation relationships (N − s > 0) in which each conservation relation represents conservation of a chemical unit or moiety the (mass) conservation relationships are given by a set of λ = N − s algebraic equations of the form: denes the locus of equilibria of the system (2.1) compatible with the reaction polyhedron xed by σ ∈ R λ >0 . Note that (2.4) contains N + δ equations [6].
Let us compute the matrix: Matrix G is square of dimensions N +δ. Taking into account that D c W = B T and D α W = 0, G reads: When we x the parameters α and k, we can write (2.4) in the form: where σ i is the mass conservation constant corresponding to the conserved moiety M i (the conservation constants corresponding to the remaining moieties in M are also xed), u is given by: 8) and the function f :  Starting from (2.7) we can compute the matrix Q (1.3) as: where Q is of dimensions (N + δ) × (N + δ + 1). By (2.11) the rank of G is equal to N + δ − 1. We can partition D σ i f as follows: where, as it can be deduced from (2.4), D σ i H = 0 and with ∂W j /∂σ i being -1 for j = i and zero otherwise. Vector D σ i W has thus one and only one entry dierent from zero.
Note that G can be partitioned as: and so the matrix Q reads: For any i = 1, . . . , λ, the vector D σ i f is linearly independent of the columns of D α f .
Taking into account that by (2.10) D c H is full rank, the vector D σ i f (for any i = 1, . . . , λ) is also linearly independent of the columns of D c f , and thus, since rank(G) = N + δ − 1, we have that rank(Q) = N + δ.
Therefore, according to Lemma 1, for any M i ∈ M , f (u, σ i ) = 0 denes a locally unique and smooth curve passing through (u * T , σ * i ) T .
Let v ∈ R N +δ+1 be the tangent vector to the curve at (u * T , σ * i ) T and decompose the product Qv into: where ) and g = D σ i f . Computed in (2.13), D σ i f has a strictly negative entry (corresponding to moiety M i ). Therefore, gv N +δ+1 = 0 if and only if v N +δ+1 = 0.
Remark 1 For convenience, we have proven (c * T , σ * i ) T to be a limit point for (2.7) with respect the mass conservation constant σ i , which does not preclude (c * T , β * ) T from being a limit point for (2.7) with respect to other parameter β chosen from the set of kinetic constants.
Remark 2 With a slight abuse of notation, if (c * T , σ * i ) T is a limit point for (2.7) (with k = k * ) with respect the mass conservation constant σ i (or with respect to any other parameter), we say that the system (2.4) (and the associated reaction network) has a saddle-node at c * , k * . +D=FJAH ! Sucient condition for a saddle-node in semi-diusive reaction networks Let us consider a reaction network fullling the following assumptions: A.1 Mass action kinetics.
A.2 Capacity for a positive equilibria.
The dynamics of such reaction system are given by a set of ODEs of the form:ċ = K + S to v to (c, k) (3.1) as described in the main text, where c ∈ R N is the vector of states, K ∈ R N ≥0 is the constant inow term, S to is the N × R to stoichiometric matrix of true and outow reactions and v to : Let us remind here that we refer by key species to the species that are present in the inlet. The system (3.1) is at equilibrium if: (3. 2) The Jacobian of the system is given by: where Y r is the matrix of the source complexes of true and outow reactions. Let µ ∈ R Rto be the vector containing the uxes of the true and outow reactions. We dene: as in the main text.
Let us dene the polynomial function and its counterpart in terms of the uxes as: The relationship between uxes and concentrations is given by: where k denotes here the vector of kinetic constants of true and outow reactions.
In semi-diusive networks there is an outow reaction with rate k o i c i for each species S i ∈ S , corresponding to the degradation of species S i .
We can write (3.2) in an equivalent way as: where k o i is the degradation constant of species S i ∈ S , the remaining parameters are xed, and f :

Proposition 2 Let us consider a reaction network fullling assumptions A.1,
A.2 and A.4 with dynamics described by Eq. (3.1) and involving the set of species S . Let µ ∈ R Rto >0 represent the uxes of true and outow reactions and J(µ) and p(µ) be as dened by (3.3) and (3.5), respectively.
If there is a vector µ * ∈ R Rto >0 such that: p i (µ * ) = 0 for i not being a key species, (3.10) then, there are strictly positive c * and k * satisfying µ * = v to (c * , k * ), and a species Provided that µ * ∈ R Rto >0 fullls (3.8-3.10), we have that: and therefore c * and k * together with K * dene an equilibrium point of (3.1).
Let us take a species S i ∈ S , and express Eq.
where the rank of J evaluated at the steady state given by c * and k * is N −1, according to (3.11).
We have: where ∂f j /∂k o i is strictly negative (equal to −x i ) for j = i and zero otherwise, i.e., the vector D ko i f evaluated at a strictly positive equilibrium has one and only one entry dierent from zero (corresponding to species S i ).
Let us choose a species S i ∈ S such that rank(Q) = N at (c * T , k o i * ) T . Note that the N vectors D ko i f for i = 1, . . . , N evaluated at the strictly positive equilibrium given by c * and k * span R N and therefore, for semi-diusive networks we can always choose a species S i ∈ S such that for Let v ∈ R N +1 be the tangent vector to the curve at (c * T , k * o i ) T and decompose the product Qv into: ) and g = D ko i f . The vector g, computed in Eq. (3.15), has one and only one entry which is dierent from zero (the strictly negative entry corresponding to species S i ). Therefore, gv N +1 = 0 if and only if v N +1 = 0.
On the other hand, according to Lemma 2, Qv = 0, with v ∈ R N +1 being unique (modulo scaling).
The rank of J is N − 1 by (3.11), and therefore a vector in its kernel is also unique (modulo scaling). Then, Consequently, and since we also know that gv N +1 = 0 if and only if v N +1 = 0, we can deduce that v N +1 = 0. Therefore, according to Denition 2, Remark 3 For convenience, we have proven (c * T , k * o i ) T to be a limit point for (3.7) with respect to the degradation constant k o i , which does not preclude (c * T , β * ) T from being a limit point for (3.7) with respect to any other parameter β chosen from the set of kinetic constants (including inow, true and the remaining outow reactions).
Remark 4 With a slight abuse of notation, if (c * T , k * o i ) T is a limit point for (3.7) (with k = k * ) with respect the degradation constant k o i (or with respect to any other parameter), we say that the system (3.7) (and the associated reaction network) has a saddle-node at c * , k * .

+D=FJAH "
Interferon-receptor complex The complexes are labeled such that C 1 = R1I, C 2 = R2I, C 3 = R1R2I, The molecularity matrix is: and taking mass action kinetics, the vector of mass action monomials reads: The stoichiometric matrix is: and has rank(S)=3. Therefore, the dimension of the stoichiometric subspace is s = 3, and the formula for the deciency of the network gives: On the other hand, the dimension of the equilibrium manifold is: The network is, therefore, under-dimensioned (λ > δ). There are three mass conservation laws, dened by Note that the strictly positive vector ϑ = (1, 1, 1, 2, 2, 3) fullls ϑS = 0, as it corresponds to a closed network. We build the matrix Λ in which the entry (i, j) corresponds to the node C i in the linkage class L j , such that: and we get: Starting from the stoichiometric matrix Y and the matrix Λ we compute a basis for the deciency subspace using Eq. (8) in the main text, obtaining: Therefore we can write Aψ as in Eq. (9) in the main text, obtaining a set of M − ℓ linearly independent equations of the form: After substituting the elements of vector ψ by the corresponding monomials we obtain the following expression for the locus of equilibria: Starting from these equations we obtain: From Eq. (13) in the main text we get the following matrix G: . We formulate an optimization problem to nd parameter vectors causing the determinant of G to vanish. The objective function is F def = det(G) 2 and we use the following decision vector: We solve the optimization problem (Eq. 14 in the main text) subject to positive concentration vectors (and the equilibrium manifold equations) using the global scatter search algorithm by Egea et al. [7]. We use the implementation of the enhanced scatter search solver (eSS) in the MEIGO toolbox [8], which is freely available at http://www.iim.csic.es/ gingproc/meigo.html.
No parameter vector was found for which the objective function becomes zero. The Matlab code is provided in S1 le (Matlab_code/IFN_receptor/Variant1).

Ternary complex formation network with IFN in excess
There are N = 5 species participating in the network with concentrations and the vector of mass action monomials: The stoichiometric matrix is: The rank of the matrix S (and, therefore, the dimension of the stoichiometric subspace) is s = 3. We compute the deciency according to the formula given in the main text as: The dimension of the equilibrium manifold is: The network is, as in the previous case, under-dimensioned. There are two mass conservation laws, with matrix B given by: Starting from the stoichiometric matrix Y and the matrix Λ we compute a basis for the deciency subspace as: and, therefore, we can express Aψ as: Substituting the elements of vector ψ by the concentration monomials we obtain the equations dening the locus of equilibria: Then we can express the concentrations as: From Eq. (12) in the main text we get the matrix G: We formulate an optimization problem to nd parameter vectors for which the determinant of G vanishes. The objective function is F def = det(G) 2 and we use the following decision vector: x = (k 41 , k 14 , k 25 , k 52 , k 63 , k 36 , k 37 , k 73 , α 1 , c 4 ).
We solve the optimization problem (Eq. 14 in the main text) subject to positive concentration vectors (and the equilibrium manifold equations) using the global scatter search solver by Egea et al. [7]. No parameter vector was found for which the objective function becomes zero. The Matlab code is provided in S1 le (Matlab_code/IFN_receptor/Variant2).

Semi-diusive network with constant IFN inow
The system in Fig. 3B, taking into account all the complexes and uxes, including the inow and outow uxes indicated by dashed arrows, is a semidiusive network. The molecularity matrix Y containing the molecularity of each species in each complex is given by Eq. (4.1). The matrix containing the molecularities of the species in the source complexes for true and outow reactions of the semi-diusive network is: , and the corresponding stoichiometric matrix can be computed starting from the stoichiometric matrix of the true reaction subnetwork (closed system in Section 1.1) as: There are in total 14 uxes corresponding to true and outow reactions (we denote the ux vector by µ). From Eq. (17) in the main text we have that: and in terms of the uxes: The vector containing the basal formation rates of the key species (I, R1 and R2) is: where v 15 = k 15 , v 16 = k 16 and v 17 = k 17 are the (constant) reaction rates of the reactions r 15 , r 16 and r 17 .
Let us build a matrix Y r with columns corresponding to the source complexes for the true and outlet reactions, such as: Now, we are ready to dene the objective function F inj = S to diag(µ)Y T r to minimize, searching values of the decision vector of uxes µ that make it vanish. The search is subject, on the one hand, to the following inequality constraints (corresponding to the key species I, R1, R2): where the expressions at the left hand side of the inequalities are given by p (1), p(2) and p(3) in (4.3) and, on the other hand, to the following equality constraints (corresponding to the remaining species): where the expressions at the left hand side are given by p(4), p (5) and p (6) in (4.3). We solve the optimization problem (Eq. 18 in the main text) using the global scatter search solver by Egea et al. [7], and no positive vector of uxes was found for which the objective function becomes zero. The Matlab code is provided in S1 le (Matlab_code/IFN_receptor/Variant3).

Semi-diusive network with IFN in excess
The network in Fig. 3C, taking into account all the complexes and uxes (including the inow and outow uxes indicated by dashed arrows) is a semi-diusive network. The molecularity matrix Y containing the molecularity of each species in each complex is given by Eq. (4.2). The matrix containing the molecularities of the species in the source complexes for true and outow reactions of the semi-diusive network is: , and the stoichiometric matrix can be computed starting from the molecularity matrix of its true-reaction network counterpart (closed network in Section 1.2): There are in total 13 uxes corresponding to true and outow reactions (contained in the vector µ). From Eq. (17) in the main text we have that: and in terms of the uxes: Note that the equality constraints can be used to express some uxes in terms others and reduce the number of decision variables.
The vector of basal formation rates of the key species, R1 and R2, is: where v 14 = k 14 and v 15 = k 15 are the (constant) reaction rates of the reactions r 14 and r 15 . We build a matrix Y r with columns corresponding to the source complexes of the true and outlet reactions: We dene now the objective function F inj = S to diag(µ)Y T r and search for values of the decision vector of uxes µ that make it vanish. The search is subject, on the one hand, to the following inequality constraints (corresponding to the key species R1 and R2): where the expressions at the left hand side are given by p(1) and p(2) in (4.4) and, on the other hand, to the following equality constraints: where the expressions at the left hand side are given by p(3), p(4) and p (5) in (4.4). We solve the optimization problem (Eq. 18 in the main text) using the global scatter search solver by Egea et al. [7]. No positive vector of uxes was found for which the objective function becomes zero. The Matlab code is provided in S1 le (Matlab_code/IFN_receptor/Variant4).
The molecularity matrix is: The stoichiometric matrix is: The rank of the matrix S (and, therefore, the dimension of the stoichiometric subspace) is s = 6. The deciency of the network is: and the dimension of the equilibrium manifold: Starting from the stoichiometric matrix Y and the matrix Λ we compute a basis for the deciency subspace as: and, therefore, we can express Aψ as: Substituting the elements of vector ψ by the concentration monomials we obtain the equations dening the locus of equilibria: We can express the concentrations as: 3 6 c 9 = α 2 /k 4 12 c 2 = (−2α 1 /k 10 2 + α 2 /k 10 2 )/c 6 The matrix G, from Eq. (13) in the main text, reads: We have that N + δ = 11, which gives us the dimensions of the matrix G (11 × 11). We formulate the optimization problem dened by Eq. (14) in the main text, to nd parameter vectors for which the determinant of matrix G becomes zero. We minimize the objective function F def = det(G) 2 using as decision vector: x = (k 71 , k 17 , k 18 , k 10 2 , k 2 9 , k 11 3 , k 3 6 , k 5 4 , k 4 12 , α 1 , α 2 , c 1 ).
We solve the optimization problem subject to positive concentration vectors (and the equilibrium manifold equations) using the global scatter search solver by Egea et al. [7]. We found a decision vector for which the objective function vanishes (see Table 5.1). At this point, the system undergoes a saddle-node bifurcation. We start a bifurcation analysis from this steady state, conrming that in fact the network has capacity for bistability (see Fig  5.1). The Matlab code is provided in S1 le (Matlab_code/STATS_early/Variant1).

Semi-diusive network
We consider the network in Fig. 5A taking into account all the uxes (including the inow and outow uxes represented by dashed arrows). The molecularity matrix is given by Eq. (5.1). We build the stoichiometric matrix of the true and outow reactions starting from matrix S for the true reaction subnetwork (closed network in section 2.1) as: There are in total 18 uxes corresponding to true and outow reactions (contained in the vector µ).
From Eq. (17) in the main text we have that: and in terms of the uxes: The vector of basal formation rates of the key species (R * , S1 and S2) is: where v 19 = k 19 , v 20 = k 20 and v 21 = k 21 are the (constant) reaction rates of the reactions r 19 to r 21 . We build a matrix Y r with columns corresponding to the source complexes of the true and outlet reactions: We can now dene the objective function F inj = S to diag(µ)Y T r , and search for values of the decision vector of uxes µ that make it vanish. The search is subject, on the one hand, to the following inequality constraints (corresponding to the key species): where the expressions at the left hand side are given by p(1), p(2) and p (3) in Eq. (5.2) and, on the other hand, to the following equality constraints: where the expressions at the left hand side are given by p(4) to p (6) in Eq. (5.2). We solve the optimization problem (Eq. 18 in the main text) using the global scatter search solver by Egea et al. [7]. We found a decision vector for which the objective function vanishes (see Table 5.2). We now compute a steady state concentration vector and a set of kinetic parameters compatible with these uxes. We x the values of the concentrations such that c = 1, and the kinetic rate constants remain equal to the uxes, i.e., k i = µ i for i = 1, . . . , 18. Starting a continuation of equilibria from this point, we conrm that the system undergoes a saddle-node bifurcation, and that it has capacity for bistability (see Fig 5.2). The Matlab code is provided in S1 le (Matlab_code/STATS_early/Variant2). Note that the equality constraints can be used to express some uxes in terms of others reducing the number of decision variables.
From Eq. (17) in the main text we have that: and in terms of the uxes: The vector of basal formation rates of the key species, R * , S1, S2, IRF 1 and CBP is: , v 28 = k 28 and v 29 = k 29 are the (constant) reaction rates of the reactions r 25 to r 29 . We build a matrix Y r with columns corresponding to the source complexes of the true and outlet reactions: We can dene now the objective function F inj = S to diag(µ)Y T r , and search for values of the decision vector of uxes µ that make it vanish. The search is subject, on the one hand, to the following inequality constraints (corresponding to the key species): where the expressions at the left hand side are obtained from p(1), p(2), p(3), p(10) and p(11) and, on the other hand, to the equality constraints: where the expressions at the left hand side are obtained from p(4), p(5), p (6), p(7), p (8), p(9) and p(12). We solve the optimization problem (Eq. 18 in the main text) using the global scatter search solver by Egea et al. [7]. We found a decision vector for which the objective function vanishes (see Table 6.1). We now compute a steady state concentration vector and a set of kinetic parameters compatible with these uxes. We x the values of the concentrations such that c = 1, and the kinetic rate constants remain equal to the uxes, i.e., k i = µ i for i = 1, . . . , 24. Starting a continuation of equilibria from this point, we conrm that the system undergoes a saddlenode bifurcation, and that it has capacity for bistability (see Fig 6.1).
We found a decision vector for which the objective function vanishes (see Table 7.1). At this point, we start a continuation of equilibrium and conrm that the system undergoes a saddle-node bifurcation. In fact, the network has capacity for bistability (see The Matlab code is provided in S1 le (Matlab_code/IFN_merged). The IFN-STAT signaling network model which we use for the analysis is compatible with the STAT activation dynamics observed upon IFN stimulation, as it can be deduced from the t of the model to experimental data ( Fig. 8.1). Experimental data were provided by Ignacio Moraga and Sandra Pellegrini from the IFNaction Consortium.