A method for inverse bifurcation of biochemical switches: inferring parameters from dose response curves

Background Within cells, stimuli are transduced into cell responses by complex networks of biochemical reactions. In many cell decision processes the underlying networks behave as bistable switches, converting graded stimuli or inputs into all or none cell responses. Observing how systems respond to different perturbations, insight can be gained into the underlying molecular mechanisms by developing mathematical models. Emergent properties of systems, like bistability, can be exploited to this purpose. One of the main challenges in modeling intracellular processes, from signaling pathways to gene regulatory networks, is to deal with high structural and parametric uncertainty, due to the complexity of the systems and the difficulty to obtain experimental measurements. Formal methods that exploit structural properties of networks for parameter estimation can help to overcome these problems. Results We here propose a novel method to infer the kinetic parameters of bistable biochemical network models. Bistable systems typically show hysteretic dose response curves, in which the so called bifurcation points can be located experimentally. We exploit the fact that, at the bifurcation points, a condition for multistationarity derived in the context of the Chemical Reaction Network Theory must be fulfilled. Chemical Reaction Network Theory has attracted attention from the (systems) biology community since it connects the structure of biochemical reaction networks to qualitative properties of the corresponding model of ordinary differential equations. The inverse bifurcation method developed here allows determining the parameters that produce the expected behavior of the dose response curves and, in particular, the observed location of the bifurcation points given by experimental data. Conclusions Our inverse bifurcation method exploits inherent structural properties of bistable switches in order to estimate kinetic parameters of bistable biochemical networks, opening a promising route for developments in Chemical Reaction Network Theory towards kinetic model identification. Electronic supplementary material The online version of this article (doi:10.1186/s12918-014-0114-2) contains supplementary material, which is available to authorized users.


.1 Network graph and ODE system
In this section we include the derivation of the set of ODEs for the protein activation network starting from the graph of complexes. The mechanism: describes the activation of the protein S mediated by the protein E where S represses its own activation. Consider also the degradation of the proteins S and P , and the constitutive formation of S, described by the following reactions: . The number of complexes is n = 8 and the graph of complexes is included in the main text ( Fig. 1), where every complex has been labeled with a number. Note that the labeling of the complexes is arbitrary. The graph consists of ℓ = 3 linkage classes. In this case we chose: C 1 = ES, C 2 = ESS, C 3 = S, C 4 = E + S, C 5 = P , C 6 = ES + S, C 7 = ∅, C 8 = E + P . The molecularity matrix Y ∈ R m×n as dened in the main text, contains the molecularities of the species involved in every complex: (1.1) The vector ψ(c) ∈ R n ≥0 containing the mass action monomials associated to the complexes (Eq. 1 in the main text) is: 2) and the reaction rates of the network are shown in Table 1.1.
The matrix A ∈ R n×n is: and the set of ODEs describing the dynamics (corresponding to Eq. 2 in the main text) reads:

Mass conservation laws and DAE system
In this section we include the derivation of the DAE system for the protein activation network. The stoichiometric subspace is spanned by the following vectors: The left null space of Y A is: The function W (dened by Eq. 3 in the main text) reads: and the reaction polyhedron consists of all c ≥ 0 such that W (c; c 0 ) = 0. In other words, there is one conservation law in the system (Eq. 4 in the main text) given by: Note that the right hand side of this expression is the total amount of enzyme E T : The existence of conservation laws implies that the ODE system (1.3) is not minimal. We can partition the state space to get a minimal set of ODEs and a set of algebraic equations, obtaining the following DAE system consisting of 4 linearly independent ODEs and one algebraic equation:

Deciency subspace and tangent bifurcation condition
In this section we compute the deciency subspace, the equilibrium manifold equations and the tangent bifurcation condition for the protein activation network. The deciency of the network is one and, therefore, the deciency subspace is one dimensional. In order to compute a basis for the deciency subspace let us rst dene a n × ℓ matrix Λ where the entry (i, j) corresponds to the node C i in the linkage class L j , being: as indicated in the main text. The matrix Λ results: A basis ω is computed taking into account Eq. 6 in the main text, obtaining: At equilibrium the following equality (Eq. 7 in the main text, with δ = 1) must hold: obtaining: Using (1.2) we can write the previous set of equations in terms of the concentrations: (1.12) The locus of equilibria in the space of the species (Eq. 9 in the main text) is given by: Let us compute the Jacobian of H s with respect to c: and the Jacobian of H s with respect to α is: (1.14) 5 The Jacobian of W with respect to c is: and the Jacobian of W with respect to α is 0. From Eq. 10 in the main text, we get: 16) and the matrix G, given by Eq. 11 from the main text reads: Computing the determinant of the matrix G, Eq. 12 from the main text (condition for multistationarity) reads: +D=FJAH Gene regulation toggle switch

Network graph and ODE system
In this section we include the derivation of the set of ODEs for the gene regulation switch example starting from the graph of complexes. The gene regulatory network under study consists of the following set of reactions: In total, there are seven species (m = 7) involved in ten reactions (r = 10). Let us order the species such that the concentrations are The number of complexes is n = 13, and the graph of complexes is depicted in Fig. 7 from the main text, where every complex has been labeled with a number. Note that the labeling of the complexes is arbitrary. The graph consists of ℓ = 6 linkage classes. In this case we select the labels indicated in Table 2.1. The molecularity matrix Y ∈ R 7×13 contains the molecularities of the species involved in every complex: (2.1) The vector ψ(c) ∈ R n ≥0 , as dened by Eq. 1 in the main text, contains the mass  v 3 10 k 3 10 c 3 r 3 v 9 4 k 9 4 c 1 c 7 r 4 v 4 9 k 4 9 c 4 r 5 v 11 5 k 11 5 c 3 c 2 r 6 v 5 11 k 5 11 c 5 r 7 v 12 6 k 12 6 c 5 c 2 r 8 v 6 12 k 6 12 c 6 r 9 v 2 13 k 2 13 c 2 r 10 v 7 13 k 7 13 c 7 action monomials associated to the complexes: ψ(c) = (c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , c 7 , c 1 c 2 , c 1 c 7 , c 3 c 7 , c 2 c 3 , c 2 c 5 , 1) T , (2.2) and the reaction rates of the network are shown in Table 2.2. The matrix A ∈ R n×n is: , and the set of ODEs corresponding to Eq. 2 in the main text reads: c 1 = k 4 9 c 4 − k 9 4 c 1 c 7 c 2 = k 1 8 c 1 − k 2 13 c 2 + k 5 11 c 5 + k 6 12 c 6 − k 11 5 c 2 c 3 − k 12 6 c 2 c 5 c 3 = k 5 11 c 5 − k 11 5 c 2 c 3 c 4 = −k 4 9 c 4 + k 9 4 c 1 c 7 c 5 = −k 5 11 c 5 + k 6 12 c 6 + k 11 5 c 2 c 3 − k 12 6 c 2 c 5 c 6 = −k 6 12 c 6 + k 12 6 c 2 c 5 c 7 = k 3 10 c 3 + k 4 9 c 4 − k 7 13 c 7 − k 9 4 c 1 c 7 . (2.3)

Mass conservation laws and DAE system
In this section we include the derivation the DAE system for the gene regulation switch example starting from the graph of complexes. The stoichiometric subspace is spanned by the following vectors: The dimension of the stoichiometric subspace is s = 5 and thus we have two conservation relation (m − s = 2). The left null space of Y A is: The existence of conservation laws implies that the ODE system given by (2.3) is not minimal. We can partition the state space to get a minimal set of 5 linearly independent ODEs and a set of 2 algebraic equations obtaining the following DAE system: c 1 = k 4 9 c 4 − k 9 4 c 1 c 7 c 2 = k 1 8 c 1 − k 2 13 c 2 + k 5 11 c 5 + k 6 12 c 6 − k 11 5 c 2 c 3 − k 12 6 c 2 c 5 c 3 = k 5 11 c 5 − k 11 5 c 2 c 3 c 5 = −k 5 11 c 5 + k 6 12 c 6 + k 11 5 c 2 c 3 − k 12 6 c 2 c 5 c 7 = k 3 10 c 3 + k 4 9 c 4 − k 7 13 c 7 − k 9 4 c 1 c 7 .

Deciency subspace and tangent bifurcation condition
In this section we compute the deciency subspace, the equilibrium manifold equations and the tangent bifurcation condition for the gene regulatory network. The deciency of the network is two δ = 13 − 6 − 5 = 2 and, therefore, the deciency subspace is two dimensional. In order to compute a basis for the deciency subspace let us rst build the n × ℓ matrix Λ where the entry (i, j) corresponds to the node C i in the linkage class L j , being as indicated in the main text. The matrix Λ results: A basis ω is computed taking into account Eq. 6 in the main text, obtaining: , (2.14) and the Jacobian of W with respect to α is the zero vector D α W = (0, 0). Therefore, the matrix G given by Eq. 11 in the main text reads: (2.15) Computing the determinant of the matrix G and applying Eq. 12 in the main text, we obtain the condition for multistationarity.