Sign-sensitivities for reaction networks: an algebraic approach

This paper presents an algebraic framework to study sign-sensitivities for reaction networks modeled by means of systems of ordinary differential equations. Specifically, we study the sign of the derivative of the concentrations of the species in the network at steady state with respect to a small perturbation on the parameter vector. We provide a closed formula for the derivatives that accommodates common perturbations, and illustrate its form with numerous examples. We argue that, mathematically, the study of the response to the system with respect to changes in total amounts is not well posed, and that one should rather consider perturbations with respect to the initial conditions. We find a sign-based criterion to determine, without computing the sensitivities, whether the sign depends on the steady state and parameters of the system. This is based on earlier results of so-called injective networks. Finally, we address systems with multiple steady states and the restriction to stable steady states.


Introduction
One of the main challenges in molecular and systems biology is to infer mechanistic details of the processes that underlie available experimental data. A common strategy consists in applying controlled perturbations to the system of interest, and compare model predictions with gathered quantitative data. As an example, consider a simple two-component system in which a histidine kinase E transfers a phosphate group to a response-regulator S: Let us imagine a protein I, which forms an inhibitory complex Y with the substrate by binding exactly one of the two forms, S or S p . This gives rise to two rivaling inhibition models with the following additional reactions In this small (and artificial) example, the measurement of the concentration of I at steady state for two starting conditions that only differ slightly in the amount of kinase E, helps us to qualitatively discriminate between these two models. Indeed, the derivative of the concentration of I at steady state increases with E in the first model, while it decreases in the second model. This exemplifies the basis of perturbation-based studies, in which the response of a system to an intervention, that being the addition of a protein, knockdown of a component, modification of reaction rate constants etc, is recorded [12,16,23,30]. In this paper we focus on how to predict this response given the model, such that comparison with experimental data can be in place.
The modeling setting is based on reaction networks and their associated evolution equations for the concentrations of the species in the network. In the examples, we employ the mass-action assumption, although this is not required for the theoretical framework. We consider perturbations on the parameters of the model, these typically being either initial concentrations or total amounts, and kinetic parameters.
We assume that the system is at steady state, and that a small perturbation to the parameter vector is performed. If the perturbation is small enough and the steady state is non-degenerate and stable, then the system is expected to converge to a new steady state. Our goal is to determine, component-wise, the sign of the difference of the two steady states as a function of the starting steady state and the parameters of the system. Mathematically, this translates into determining the signs of the derivatives of all concentrations at steady state with respect to the performed perturbation; these signs are referred to as sign-sensitivities [27].
Numerous computer-based approaches exist to find sign-sensitivities, with the same or different modeling framework as the one we use here [15,16,29]. A general strategy assumes all parameters known except the one of interest, and the response of the system is investigated using simulations. However, most quantities and parameters are notoriously hard to measure or estimate, which reduces the applicability of these methods. It is worth highlighting the algorithm by Sontag [27] to determine sign-sensitivities with respect to the increase of a total amount (a quantity that is preserved under the dynamics of the system), while treating the rest of the parameters as unknowns. Alternatively, if the network is small enough, one can attempt direct manipulation of the steady state equations [8,11]. In recent works [20,21], Okada and Mochizuki provide a theorem to determine zero sensitivities from network structure alone. Similarly, Brehm and Fiedler study whether the sensitivity is zero [2]. Another series of works [24,25] address the question of "how large the absolute value of the sensitivity is", by finding upper bounds for reaction networks in a specific class.
In this paper we derive a closed formula for the sensitivities. When the kinetics are polynomial, then the derivative is expressed as a rational function in the parameters of the system and the concentrations of the species at the steady state. If the numerator and the denominator of this function have all coefficients with the same sign, then the sign of the derivative is easily determined, and it does not depend on the initial steady state nor on the parameters of the system. When the signs of the coefficients differ, then one can employ standard techniques, such as those based on the Newton polytope, e.g. [3], to determine whether the derivative can both be positive and negative, thereby concluding that the sign depends on the chosen steady state and/or parameter values. For example, for the two rivaling inhibitory models above, we find that the derivative c ′ I of the concentration of I with respect to the initial concentration of E at a steady state c for model 1 and model 2 are where k i > 0 stands for the reaction rate constants of the reactions in the network (in the order given above), c Z denotes the concentration of the species Z at the steady state, and By inspecting the signs of these polynomials at positive values of c and k, we conclude that model 1 leads to an increase of the concentration of I while model 2 to a decrease.
In [23], apparently paradoxical results on sign-sensitivities were brought up to attention. We recover these phenomena in this work, and encounter a new surprising counter-intuitive result: the concentration of a species X at steady state might decrease as a function of X itself; that is, the concentration of X might decrease after the addition of X to the system (Section 6).
After exemplifying how to find the sensitivities for different types of perturbations (Section 3), we focus on perturbations of concentrations. We argue that, mathematically, it is better posed to discuss responses to a change of an initial concentration rather than to a change of a total amount (Section 4). We then employ recent results relating the sign of the determinant of the Jacobian of a function with sign-vectors conditions from [18] to determine, without computing the sensitivities, whether the sign depends on the steady state and parameters of the system (Section 5). We conclude by discussing the existence of multiple steady states and the restriction to stable steady states, which are the only ones observable in an experimental setting (Section 6).

Reaction networks
The processes we consider are modeled by reaction networks, which can be seen as directed graphs. Specifically, a reaction network consists of a set of species {X 1 , . . . , X n }, and a directed graph whose nodes are finite linear combinations of species (called complexes). The directed edges are called reactions. We let r be the number of reactions. An example of a reaction network [27], modeling the transfer of phosphate groups from a kinase E to a substrate S that has two phosphorylation sites is: The species of the network are E, E p , S 0 , S 1 , S 2 : E, E p are the unphosphorylated and phosphorylated forms of the kinase E and S 0 , S 1 , S 2 denote the substrate with no, one or two phosphate groups attached. The complexes are S 0 + E p , S 1 + E, S 1 + E p and S 2 + E. This network, which is used as a running example, originates from the work of Sontag on sign-sensitivities as well [27]. See [7,14] for an expanded introduction to reaction networks. The source of a reaction is called the reactant, while the target is called the product. We assume the set of species is numbered such that each complex y can be identified with a vector in R n ; for instance, the complex X 1 + 2X 2 is identified with the vector (1, 2, 0, . . . , 0) in R n , where n is the number of species. In this way, each reaction y → y ′ gives rise to a vector y ′ − y in R n , encoding the net production of each species with respect to the reaction. After choosing an order of the set of reactions, these vectors are gathered as columns of a matrix, called the stoichiometric matrix N ∈ R n×r . The stoichiometric matrix for network (1) is where the species are ordered as S 0 , S 1 , S 2 , E, E p .
As it is custom within chemical reaction network theory, in this work we model the evolution of the concentration of the species in the network in time by means of a system of ordinary differential equations (ODEs). Specifically, denote by x i (t) the concentration of X i at time t (or x C (t) for a species C). One chooses the rate v y→y ′ (x) of each reaction of the network, to be a differentiable function from R n ≥0 to R ≥0 , and gathers these rates into a vector v(x) from R n ≥0 to R r ≥0 , using the established orders of the sets of species and reactions. Then the evolution equations of the vector of concentrations x = (x 1 , . . . , x n ) takes the form The vector v(x) depends often on parameters, as it is exemplified below. Therefore, we often write v k (x) to indicate dependence on some vector of parameters k, and define Under the assumption of mass-action kinetics, we have x yi i , with 0 0 = 1. Here, k y→y ′ > 0 is called the reaction rate constant and is treated as a parameter, since it is often unknown. Let B be the n × r matrix such that the i-th column is the reactant vector of the i-th reaction. Then, under the assumption of mass-action kinetics, the right-hand side of system (3) can be written equivalently as The reaction rate constant is incorporated in the reaction network as a label of the reactions, such that we write for the network in (5) We write x 1 , . . . , x 5 for the concentrations of S 0 , S 1 , S 2 , E, E p respectively. Under the assumption of mass-action kinetics, the matrix B and the vector v(x) are which together with the matrix N in (2), give the following ODE system: It is clear from (3) that the vector dx dt belongs to the column span Im(N ) of N , called the stoichiometric subspace. Thus, given an initial condition x 0 , the solution to (3) is confined to the linear subspace x 0 + Im(N ). Further, both R n ≥0 and R n >0 are also forward-invariant for the trajectories of (3) [26]. Each of the sets (x 0 + Im(N )) ∩ R n ≥0 ⊆ R n ≥0 is called a stoichiometric compatibility class. In this work we parametrize these sets in two ways. First, choose a matrix W whose rows form a basis of Im(N ) ⊥ . Then, the set (x 0 + Im(N )) ∩ R n ≥0 agrees with the subset of R n ≥0 defined by the equation x ∈ R n ≥0 . This set is independent of the choice of matrix W and is parametrized by x 0 ∈ R n ≥0 . Let now d be the dimension of Im(N ) ⊥ . Alternatively, one might consider vectors T = (T 1 , . . . , T d ) ∈ R d and consider the sets {x ∈ R n ≥0 | W x = T }. Each such set corresponds to a stoichiometric compatibility class. However, the set depends both on W and T , that is, the vector T alone does not characterize the class. We refer to T as the vector of total amounts and to W as a matrix of conservation laws.
A matrix of conservation laws W for network (1) is which gives rise to the following equations for the stoichiometric compatibility class with vector of total amounts (T S , T p , T E ) The first equation encodes that the substrate is conserved, the third that the kinase is conserved, and the second that the phosphate group is either in S 1 , S 2 or E p , with S 2 having two sites. Alternatively, we can write the stoichiometric compatibility class of x 0 as We illustrate that it is more meaningful to use this second parametrization when studying sign-sensitivities in Section 4.
In this work, we are interested in the positive steady states of the system (3) restricted to a stoichiometric compatibility class. These are defined by a system of equations Since d equations of f k (x) are redundant (as W N = 0), we remove them from the system f k (x) = 0 and obtain a system with n equations and n variables, which we write as after removing redundant equations, and the last d components are W (x − x 0 ) tr = 0. Here tr stands for the transpose of a vector or matrix. If we wish to parametrize stoichiometric compatibility classes with T ∈ R d , then we simply write F k,T (x) = 0 for the corresponding system. For network (1), we remove the equations corresponding to the species S 0 , S 1 , E, and obtain the following system defining the steady states in the stoichiometric compatibility class of x 0 : We conclude this section with a definition: we say that a steady state x * is degenerate if the Jacobian of F k,x * (x) evaluated at x * is singular, that is, has vanishing determinant. This is equivalent to the Jacobian of the function f k (x) be singular on Im(N ), c.f. [32, Eq (6.1)].

Sign-sensitivities
We consider a reaction network with associated ODE system (3) and a steady state c. In this section we investigate how the steady state c changes upon a small perturbation to a parameter of the system, that being either k or x 0 (or T ). Specifically, we consider the vector of parameters k ∈ R m of the rate function v k (x), and the vector of parameters of initial conditions x 0 ∈ R n , or the vector of parameters of total amounts T ∈ R d . These live in a subspace Ω of R M with M = m + n or M = m + d, and we write generically α ∈ Ω for the vector of parameters of either form (k, x 0 ) or (k, T ). We let γ 0 ∈ Ω be the vector of parameters corresponding to our steady state c, such that F γ0 (c) = 0.
A formula for sensitivities. We consider a continuously differentiable map where ǫ > 0 and such that γ(0) = γ 0 . If c is not degenerate, then the Implicit Function Theorem implies that locally around 0, there is a continuously differentiable curve c(s) with c(0) = c and such that c(s) is a steady state of the reaction network and stoichiometric compatibility class with parameters γ(s), that is, F γ(s) (c(s)) = 0.
The question we address here is how to determine the sign of the derivative of c(s) with respect to s at s = 0, which we denote by c ′ (0). We let γ ′ (s) denote the derivative of γ with respect to s.
We view F α (x) as a function in both α and x and let J α,1 (x) := ∂Fα(x) ∂x ∈ R n×n denote the Jacobian matrix of F α (x) with respect to the vector x and similarly J α,2 (x) := ∂Fα(x) ∂α ∈ R n×M denote the Jacobian matrix of F α (x) with respect to the vector α. Differentiation of F γ(s) (c(s)) = 0 with respect to s and evaluation at s = 0 gives This results in a linear system in n unknowns c ′ 1 (0), . . . , c ′ n (0) with coefficient matrix J γ0,1 (c) and independent term J γ0,2 (c) · γ ′ (0), which we can find if the steady state c and γ 0 are given. Since the steady state c is non-degenerate, the coefficient matrix has full rank n, and hence this system has a unique solution. Note that neither the coefficient matrix of the linear system J γ0,1 (c) nor J γ0,2 (c) depend on the specific perturbation γ. Further, the last d rows of J γ0,1 (c) are W . When α = (k, x 0 ), the last d rows of J γ0,2 (c) are (0 | −W ), where 0 is the zero matrix of size d × m. Similarly, when α = (k, T ), the last d rows of J γ0,2 (c) are (0 | −I d×d ). In both cases the upper n − d rows are zero in the last d entries.
Using Cramer's rule, c ′ i (0) is expressed as a fraction where the denominator is the determinant of J γ0,1 (c) and the numerator is the determinant of the matrix obtained by replacing the i-th column of J γ0,1 (c) by −J γ0,2 (c) · γ ′ (0). When the rate functions are mass-action, then c ′ i (0) becomes a rational function in the parameters and the entries of c = (c 1 , . . . , c n ).
In the general scenario, nor the steady state c nor the parameter value γ 0 are known, and therefore we aim at determining the sign of c ′ i (0) for all values of c and γ 0 , and at deciding whether this sign is independent of these values. As it has been used in several works, e.g. [3,6], the set of all positive steady states is studied by means of a parametrization ϕ: U → R n >0 , such that the image of ϕ is the set of positive steady states (see [3] for strategies to find parametrizations).
We illustrate this framework and computations with selected perturbations γ for our running example (1). First, note that due to the matrix of conservation laws in (7), the steady state equations for x 1 , x 2 and x 4 are redundant. Thus, a positive steady state is simply a point x ∈ R 5 >0 satisfying the steady state equations for x 3 and x 5 : Any solution to this system is of the form that is, the set of positive steady states is parametrized by x 1 , x 2 and x 4 , where U = R 3 >0 . If these three variables are positive, then so is c.
>0 be the vector of parameters. The function F α (x) is taken to be the left-hand side of the system in (10). This leads to the following Jacobian matrices: Perturbing k 1 .
We consider first the perturbation that maps k 1 to k 1 + s. This gives γ ′ (0) = (1, 0, 0, 0, 0, 0, 0, 0, 0) and hence J α,2 (x) · γ ′ (0) is simply the first column of J α,2 (x) in (15), which is (0, −x 1 x 5 , 0, 0, 0) tr . We solve system (11) with these data and obtain for c = (x 1 , . . . , x 5 ) that Table 1. Sign-sensitivities with respect to adding a small amount to each of the parameters. Each column gives the sign-sensitivity with respect to one parameter. τ 1 is the sign of k 1 x 1 + k 2 (x 2 − x 4 ) and τ 2 is the sign of k 1 k 4 x 2 1 − k 2 k 3 x 2 2 at the steady state (which can be zero). Gray cells are determined using the parametrization, and for the other cells the sign is determined for all x ∈ R 5 >0 . where . We readily see that c 1 and c 5 decrease, and c 2 and c 4 increase when k 1 is slightly increased. The sign of c ′ 3 (0) is not determined (yet). But we have not imposed that c is a steady state. In order to do that, we evaluate c ′ 3 (0) in the parametrization and obtain that the sign of c ′ 3 (0) at a steady state is the sign of Clearly, this expression can be positive, negative or zero, after appropriately choosing k 1 , k 2 , x 1 , x 2 , x 4 . We conclude that the sign of the change of c 3 with respect to this perturbation is not parameter and variable independent, and therefore information on the specific value of the steady state is required.
Perturbing x 0 4 . We perturb now x 0 4 (that is, the initial concentration of E) by the addition of a small amount s. We now have γ ′ (0) = (0, 0, 0, 0, 0, 0, 0, 1, 0) and hence J α,2 (x) · γ ′ (0) is the eighth column of J α,2 (x), which is (0, 0, 0, 0, −1) tr . We solve the resulting system (11) and obtain for c = (x 1 , . . . , x n ) that After evaluating the numerator of c ′ 1 (0) in the parametrization, we obtain which only attains positive values. We conclude that the concentration of S 0 at steady state increases when an infinitesimal amount of E is added to the system.
We proceed in the same way to determine c ′ i (0) after perturbing each of the reaction rate constants k j and initial concentrations x 0 j one by one by adding a small amount s. The sign-sensitivities are summarised in Table 1. A seemingly striking insight of this table is that an increase of a certain species can be paired with both an increase or decrease of another species, depending on the perturbation applied. For instance, E p increases (c ′ 5 (0) > 0) after increasing either x 0 3 or x 0 4 , while E decreases (c ′ 4 (0) < 0) for the first perturbation and increases (c ′ 4 (0) > 0) for the second. This highlights that perturbation studies need to be appropriately interpreted, as it would be wrong to conclude, out of the column for x 0 3 , that E decreases when E p increases. That is, one needs to pair the direct perturbation to the response, and not two responses to a perturbation. This "paradoxical" result has been first pointed out in [23].
General perturbations. The outlined framework accommodates all types of perturbations, not only consisting in adding a small amount to one of the parameters. We illustrate this with two perturbations: in the first we scale two reaction rate constants by s, and in the second a small amount s is added to x 0 4 and x 0 5 simultaneously.
Consider next the perturbation to α such that γ(s) = (k 1 , k 2 , k 3 , k 4 , x 0 1 , x 0 2 , x 0 3 , x 0 4 + s, x 0 5 + s). Then γ ′ (0) = (0, 0, 0, 0, 0, 0, 0, 1, 1) and J γ0,2 (c) · γ ′ (0) is the sum of the last two columns of J γ0,2 (c), namely the vector (0, 0, 0, −1, −2) tr . Then the solution c ′ i (0) to the corresponding system is the sum of c ′ i (0) for the perturbation x 0 4 → x 0 4 + s and c ′ i (0) for the perturbation x 0 5 → x 0 5 + s. By Table 1, the sign of c ′ 4 (0) and c ′ 5 (0) is +. We further obtain that the sign of c ′ 1 (0), c ′ 2 (0) and c ′ 3 (0) can be any of −, 0, +. In this example it is straightforward to decide whether the numerator of c ′ i (0) can attain any sign when the polynomial has coefficients of both signs. For larger systems, an often successful approach consists on investigating the vertices of the Newton polytope associated with the polynomial. If two of the vertices correspond to monomials with coefficients of opposite signs, then the polynomial attains all signs for positive values of the variables. This strategy has been used in numerous recent works with chemical reaction network theory, e.g. [3,5,19], and we refer the reader to [3] for an expository account.
Remark 3.1. In [2] the authors provide structural conditions to determine whether a sign-sensitivity is zero, for ODE systems arising from a subclass of rate functions that does not include mass-action. In that work, perturbations on reaction rate constants of the form k i → k i + s are considered using the corresponding equation (11). In Metabolic Control Analysis (MCA) [10], so called flux/concentration control coefficients are considered. The latter measures sensitivity similarly to here, as the derivative of the logarithm of a concentration c i with respect to the logarithm of another concentration x 0 j is taken, or what is equivalent j ci(0) . These are found using (11) as well, after adjusting the formula. Since x 0 j ci(0) is positive, this factor is redundant when considering sign-sensitivities, but in MCA, of relevance is the value of this (normalized) derivative, and not only its sign. See [13] for a gentle introduction to MCA and control coefficients.
Although we did not change the expression for the total amount T p , the sign-sensitivities changed drastically. For example, S 0 decreases when T p is increased for the first matrix of conservation laws, while it decreases for the second choice. We conclude that perturbations with respect to total amounts might not be meaningful and need to be appropriately interpreted as perturbations of the considered system.
We proceed to investigate perturbations with respect to initial concentrations, and show that in this case, the sign-sensitivities do not depend on the choice of matrix of conservation laws. For the rest of the section we let α = (k, x 0 ). Note that J α,2 (x) is independent of x 0 since F α (x) is linear in x 0 . Proof. Let W, W ′ be two matrices of conservation laws. Then there exists an invertible d × d matrix A such that W ′ = AW . Let F α (x) and F ′ α (x) be the corresponding steady state functions from (10), and denote by J, J ′ (with the appropriate subindices) their Jacobian matrices respectively. Then where I n−d is the identity matrix of size n − d (see text after (11)). Since I n−d 0 0 A is invertible, the solution to (11) for W and W ′ agree.
Having established that c ′ j (0) does not depend on the choice of matrix of conservation laws, we can easily prove a series of lemmas by appropriately selecting this matrix. First, we note that if a concentration does not take part of any conservation law, then all sensitivities with respect to changes to this concentration are zero. Proof. By hypothesis, the i-th column of any matrix of conservation laws is zero. Consequently J γ0,2 (c) · γ ′ (0) is the zero vector and the only solution to (11) is the zero vector.
In the next proposition we discuss perturbations with respect to an initial concentration that only appears in one conservation law, and how it relates to the perturbation with respect to the corresponding total amount. Proposition 4.3. Assume the matrix of conservation laws W = (w j,i ) is such that the i-th column has only one non-zero entry, that is, there exists ℓ such that w ℓ ′ ,i = 0 for ℓ ′ = ℓ and w ℓ,i = 0.
Let M j be the minor of J γ0,1 (c) obtained by removing the j-th column and the (n − d + ℓ)-th row, divided by det J γ0,1 (c). Then The statement of the proposition follows after noticing that J γ0,2 (c) · γ ′ i (0) is the vector with −w ℓ,i in the (n − d + ℓ)-th entry and zero everywhere else, and the vector J γ0,2 (c) · (γ * ℓ ) ′ (0) is −1 in the (n − d + ℓ)-th entry and zero otherwise.
In particular, it follows from Proposition 4.3 that if w ℓ,i > 0, then the perturbations with respect to x 0 i or T ℓ yield sensitivities with the same sign. As a consequence, if the matrix W is row reduced, then perturbation with respect to a slight increase of a total amount can be interpreted as the perturbation with respect to x 0 i for i the index of the first non-zero entry of the corresponding conservation law. Hence, the value of the perturbation with respect to this total amount is well defined under the restriction that the other conservation laws do not involve x i . An immediate consequence of Proposition 4.3 is that if two columns i, i ′ of W have both only one non-zero entry, at the ℓ-th row, then sensitivities with respect to x 0 i and x 0 i ′ agree up to the product of these non-zero entries. This is summarized in the following lemma.
then c ′ j (0) for the perturbation γ i sending x 0 i to x 0 i + s agrees with that for the perturbation γ i ′ sending x 0 i ′ to x 0 i ′ + s times w ℓ,i ′ /w ℓ,i . Proof. By Proposition 4.3, the numerator of c ′ j (0) for the perturbation γ i is (−1) n−d+ℓ+j w ℓ,i times the minor of J γ0,1 (c) obtained by removing the j-th column and the (n − d + ℓ)-th row, and for γ i ′ , this same minor is multiplied by (−1) n−d+ℓ+j w ℓ,i ′ .
As an example, consider the two rivaling models in the introduction. For model 1, a matrix of conservation laws is with the order of species E, E p , S 0 , S 1 , I, Y . The pairs (E, E p ) and (S 0 , S 1 ) satisfy the hypothesis of Lemma 4.4. Therefore, for any j, the sign of c ′ j (0) is the same when either E or E p are increased, and similarly, is the same when either S 0 or S 1 are increased.

Parameter-independent sign-sensitivities
Although the computation of c ′ i (0) by solving system (11) might seem straightforward, it requires the computation and analysis of two symbolic determinants. As noticed in [1,9], these computations are expensive for relatively large reaction networks, as those encountered in applications. In this section we investigate an alternative approach to decide whether the sign of c ′ i (0) does not depend on the value of the parameters nor on c = (x 1 , . . . , x n ) ∈ R n >0 , that is, without imposing that c is a steady state. We do this for perturbations of an initial concentration as in the previous section. Once it has been established that the sign is independent of the parameters and x, then it can easily be determined after arbitrarily choosing values.
The subsequent results are based on the study of injective networks and sign vectors, as presented in [18]. For that, some notation needs to be introduced. The sign-vector σ(v) of a vector v is obtained by taking the sign component-wise. For V ⊆ R n , let σ(V ) be the set of sign vectors of all elements in V , and Σ(V ) be the subset of R n containing all, possibly lower dimensional, orthants that V intersects. Consider a function in R n >0 of the form with N ∈ R n×r of rank n − d, B ∈ R r×n , k ∈ R r >0 . Let S ⊆ R n be a vector subspace of dimension d and N ′ be a submatrix of N given by n − d linearly independent rows of N (such that ker(N ) = ker(N ′ )). Define a function G k (x) with the first n−d components equal to N ′ diag(k)x B , and the last d components W x tr , for W any basis of S ⊥ . Then the determinant of the Jacobian of G k is a polynomial in k and x such that all coefficients have the same sign if and only if (16) σ(ker(N )) ∩ σ(B tr (Σ(S\{0}))) = ∅.
(see [18]). Further, if any of these conditions hold, the function G k (x) is injective on all cosets (x 0 + S) ∩ R n >0 for any choice of k. To understand how (16) arises, one first notes that the determinant of the Jacobian J G k of G k is a polynomial in k and x such that all coefficients have the same sign if and only if it never vanishes. Vanishing of det J G k means that the Jacobian of G k has non-trivial kernel, that is, there exists a non-zero vector u in ker(J g k ) which further satisfies W u = 0, i.e. u ∈ S \ {0}. Using J g k (x) = N diag(k)B tr diag( 1 x ), we have u ∈ ker(J g k (x)) if and only if ker(N ) contains diag(k)B tr diag( 1 x )u. Condition (16) arises from noticing that varying x and u ∈ S \{0} means considering all orthants that S \ {0} intersects, that is, Σ(S\{0}), and then a vector k such that diag(k)B tr diag( 1 x )u belongs to ker(N ) exists if and only if B tr diag( 1 x )u has the sign of some vector in ker(N ). For details of this construction we refer the reader to [18].
When applying (16) to a reaction network with mass-action kinetics, we consider S = Im(N ), and if (16) holds, then the reaction network is said to be injective. By verifying the sign equality (16) with N the stoichiometric matrix, S = Im(N ), and B the exponent matrix in (4), we can determine whether the sign of the determinant J α,1 (x), and hence of the denominator of c ′ i (0), is constant. We emphasize that we do not impose that x is a steady state in the computations in this section.
In order to study the numerator, we interpret it as the Jacobian of a function of the form g k (x) as above, with the same coefficient matrix N , and suitable exponent matrix B j and vector space S j . Afterwards we apply (16). The specific form of these objects is given in the next proposition.
Proposition 5.1. Assume mass-action kinetics. Consider the perturbation γ i and assume that Im(N ) ⊥ contains vectors with non-zero i-th entry. Let W be a matrix of conservation laws such that the only row where the i-th component is non-zero is the first, where it takes the value 1. Let S j ⊆ R n−1 be the kernel of the vector subspace spanned by the rows of the matrix obtained from W , by deleting the first row and the j-th column, and let B j be obtained from B by removing the j-th row.
Then the sign of the numerator of c ′ j (0), as a function of k and x, is independent of k and x if and only if σ(ker(N )) ∩ σ(B tr j (Σ(S j \{0}))) = ∅. Proof. Let x = (x 1 , . . . , x j−1 , x j+1 , . . . x n ) ∈ R n−1 be the vector x where the j-th entry is deleted. Let N ′ be a matrix formed by n − d linearly independent rows of N (such that ker(N ) = ker(N ′ )), and for ℓ = 1, . . . , r, let k ℓ = k ℓ x yj j if the reactant of the ℓ-th reaction is y. With the choice of W and the considerations before Lemma 4.4, the numerator of c ′ j (0) is, up to a constant sign, the determinant of the submatrix J ′ of J α,1 (c) obtained by removing the j-th column and the (n − d + 1)-th row. This matrix J ′ agrees with the Jacobian of the function G k ( x) in R n−1 with the first n − d entries equal to N ′ diag ( k ) x Bj and bottom d − 1 entries W ′ x tr , where W ′ is the matrix obtained from W , by deleting the first row and the j-th column.
As recalled in (16), by [18], the sign of the determinant of J ′ does not depend on k nor x, hence on k nor x, if and only if the sign condition in the statement holds.
To illustrate this result, we consider network (1), the perturbation of x 0 2 by adding s, and focus on c ′ 1 (0). By Table 1, we already know that the sign of c ′ 1 (0) is +, and this holds for any x without imposing the steady state condition. In particular, in the notation of Proposition 5.1, i = 2 and j = 1.
The kernel of N in (2) is generated by the vectors (1, 1, 0, 0) and (0, 0, 1, 1), and hence for any u = (u 1 , u 2 , u 3 , u 4 ) in ker(N ), the sign of u 1 and u 2 agree, and the sign of u 3 and u 4 agree. The matrix B 1 obtained by removing the first row of B in (6) and a matrix of conservation laws satisfying the hypothesis of Proposition 5.1 are given as  (0, a, a, 0), which has the sign of a vector in ker(N ) only if a = 0, a contradiction.
We have therefore verified that the sign condition in Proposition 5.1 holds, and therefore c ′ 1 (0) does not depend neither on k nor on x. Clearly, finding the sign vectors manually is not optimal at all. In [18], see also [22], strategies to verify whether this equality holds are presented. Figure 1. (a) A simple network of a hybrid histidine kinase, taken from [17]. (b) Signsensitivities with respect to an increase of x 0 1 and x 0 5 . ± * means that the sign is + when k 1 ≥ k 3 , that is, when the network has exactly one positive steady state.

Stable vs unstable steady states
In the previous sections we have not taken into consideration whether the steady states are asymptotically stable or unstable. In practice, in an experimental setting, only stable steady states are observable. Although it is often not possible to restrict parametrizations of the set of steady states to only stable steady states, some relevant information can be extracted from the sign of the determinant of the Jacobian J α,1 (x).
Specifically, assume the function F α (x) is constructed from a matrix of conservation laws W that is row reduced, and let i 1 , . . . , i d be the indices of the first non-zero entries of the rows of W . The first n − d entries of F α (x) can be chosen to be the entries of f k (x) with index different from i 1 , . . . , i d . Let τ = d ℓ=1 (n − ℓ − i ℓ ) be the sign of the permutation that reorders the entries of F α (x) such that the entries defined by W are at entries i 1 , . . . , i d . Then, by [32,Prop. 5.3], the determinant of J α,1 (x) is (−1) τ times the product of the n − d nonzero eigenvalues of the Jacobian of f k evaluated at x. Hence, if the steady state is hyperbolic and asymptotically stable, then necessarily the sign of this determinant is (−1) τ +n−d .
In the previous examples, the determinant of the Jacobian J α,1 (x) at a steady state had a constant sign, which was actually (−1) τ +n−d , and hence in accordance with stability. We now illustrate by means of an example what can be said about sensitivities when the network has unstable steady states.
For that, we consider a simple model of a hybrid histidine kinase from [17]. The network is depicted in Figure 1(a). Under mass-action kinetics, there exist stoichiometric compatibility classes for which this network has three positive steady states [17], two of which are asymptotically stable [28]. Further, by [3], the network admits three positive steady states in some stoichiometric compatibility class if and only if k 3 > k 1 . If k 1 ≥ k 3 , then the network has exactly one positive steady state in each stoichiometric compatibility class.
We order the species as HK, HK p0 , HK 0p , HK pp , RR and RR p , and let x 1 , . . . , x 6 denote their concentrations respectively. Following [3], the set of positive steady states admits a parametrization in terms of x 1 and x 5 obtained by solving the steady states equations of x 2 , x 3 , x 4 , x 6 in these variables: We choose the matrix of conservation laws W = 1 1 1 1 0 0 0 0 1 1 1 1 , and construct F k,x 0 (x) with first four components equal to f 2 , f 3 , f 4 , f 6 . The determinant of J α,1 (x) evaluated at the parametrization yields: Here we see that if k 1 ≥ k 3 , then this determinant has negative sign, which is actually the sign it attains when the steady state is asymptotically stable and hyperbolic. Indeed, in this case n − d = 4 and (−1) τ = (−1) n−1−1+n−2−5 = (−1) 5 = −1. If k 3 > k 1 , then the stable steady states will necessarily satisfy that the sign of det J α,1 (x 1 , x 5 ) is negative. Using this, we proceed as above to compute the sign of the sensitivities with respect to adding a small amount to each of x 0 i . By Lemma 4.4, it is enough to compute the sensitivities with respect to perturbing x 0 1 (which agrees with the perturbations with respect to x 0 2 , x 0 3 and x 0 4 ) and x 0 5 (which agrees with x 0 6 ). Figure 1(b) shows the obtained sign-sensitivities under the assumption that det J α,1 (x 1 , x 5 ) is negative. If this determinant is positive, then all signs are reversed, but this implies that the steady state is unstable.
An apparently surprising property of this network is that the addition of HK 00 , that is, x 0 1 , might lead to the decrease of HK 00 . To have a closer inspection at this phenomenon, using the parametrization, we have that c ′ 1 (0) at the steady state defined by x 1 , x 5 is c ′ 1 (0) = k 2 k 5 (k 1 k 3 x 1 − k 4 k 6 x 2 5 ) / det J α,1 (x 1 , x 5 ).
By letting k 1 = · · · = k 6 = 1, the system has exactly one positive steady state in each stoichiometric compatibility classes and det J α,1 (x 1 , x 5 ) < 0. The x 1 -component of the steady state defined by x 1 = 2, x 5 = 1 will decrease after a small amount of x 1 is added to the system. On the other hand, the x 1 -component of the steady state defined by x 1 = 1, x 5 = 2 will increase.
Two-site sequential and distributive phosphorylation cycle. We conclude with one extra example where we analyze sign-sensitivities of a classical model. We consider the reaction network in which a substrate S becomes doubly phosphorylated by a kinase E and dephosphorylated by a phosphatase F . We let S 0 , S 1 , S 2 be the three phosphoforms of S with 0, 1, 2 phosphorylated sites, respectively. This gives rise to the following reactions [4,31]: We order the species as E, F, S 0 , S 1 , S 2 , ES 0 , F S 1 , ES 1 , F S 2 and let x 1 , . . . , x 9 denote their concentrations respectively. A matrix of conservation laws is The set of positive steady states admits a positive parametrization in x 1 , x 2 , x 3 , obtained by solving the system f 4 = · · · = f 9 = 0 in x 4 , . . . , x 9 , where f i is the mass-action evolution equation for x i . It is well known that this network admits between one and three positive steady states in each stoichiometric compatibility class, e.g. [4,31].
We consider det J α,1 (x) evaluated in the parametrization, and assume it is negative: namely this is the sign of this determinant when the steady state is asymptotically stable and hyperbolic. Under this assumption, we compute the sign-sensitivities c ′ 1 (0), . . . , c ′ 5 (0) with respect to a small increase of x 0 1 , x 0 2 and x 0 3 , and obtain that none of them is given by a rational function with numerator of fix sign, indicating that all signs might be possible for this system. However, it is not straightforward to analyze the sign of the numerators while imposing that det J α,1 (x) is negative.