Uncovering bifurcation behaviors of biochemical reaction systems from network topology

TL;DR

This paper introduces a method to study how biochemical networks produce different cell states by analyzing their structure, applied to macrophage polarization.

Contribution

A new method called Structural Bifurcation Analysis (SBA) is introduced to study bifurcation behaviors from network topology alone.

Findings

SBA identifies substructures and parameters responsible for bifurcation behaviors in biochemical networks.

The method was applied to macrophage polarization, revealing how network structure constrains polarization patterns.

Gene deletions like SOCS3 are predicted to cause intermediate polarization states, offering testable hypotheses.

Abstract

The regulation of biological functions is achieved through the modulation of biochemical reaction network dynamics. The diversity of cell states and the transitions between them have been interpreted as bifurcations in these dynamics. However, due to the complexity of networks and limited knowledge of reaction kinetics, bifurcation behaviors in biological systems remain largely underexplored. To address this, we developed a mathematical method, Structural Bifurcation Analysis (SBA), which decomposes the system into substructures and determines important aspects of bifurcation behaviors—such as substructures responsible for bifurcation conditions, bifurcation-inducing parameters, and bifurcating variables—solely from network topology. We establish a direct relationship between SBA and classical bifurcation analysis, enabling the study of systems even in the presence of conserved quantities. Additionally, we provide a step-by-step bifurcation analysis for general use. We applied our method to the macrophage M1/M2 polarization system. Our analysis reveals that the network structure strongly constrains possible patterns of polarization. We also clarify the dependency of the M1/M2 balance on gene expression levels and predict the emergence of intermediate polarization patterns under gene deletions, including SOCS3, which are experimentally testable.

Figures

Summary of Structural Bifurcation Analysis (modified from^[13](#CR13)^). One can construct a matrix $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document}$ according to the network structure of a given chemical reaction networks. With the matrix, one can further conclude the conditions where bifurcations can occur, the inducing parameters, and the chemicals that exhibit the bifurcations.

Bifurcation diagrams of the example network with different ODE models. The three models respectively show a (a) pitchfork, (b) transcritical, and (c) saddle-node bifurcation. Independent of the kinetics, we verify our prediction that the bifurcation behavior is exhibited only by *B*. The ODE models and the parameter settings are given by SI Appendix B.

Summary of the standard procedure of Structural Bifurcation Analysis. The theory of SBA has been extended to networks with conserved quantities allowed, and we further improve the procedure for users to facilitate the roles and relations of subnetworks in the bifurcation dynamics.

Macrophage polarization network for the wild type. By an arrow " $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dashv$$\end{document}$ " pointing from a complex to a reaction, we mean a negative regulation. The superscript $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\o

Bifurcation analysis of the WT network. (a) The bifurcation diagrams with the conserved quantity $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1=\text {STAT1}+\text {STAT1}^+$$\end{document}$ ranging in [6, 8.5] being the bifurcation parameter. Markers in different shapes are used to denote points on distinct equilibrium branches. The short and long stable equilibrium branches can be interpreted as the M1

The influence graphs and the rearranged $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document}$ matrices corresponding to the wild, the SOCS3-deleted, and the SOCS3-/KLF4- deleted types. In the rearranged matrix, we color the block $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepac

Bifurcation analysis of the modified networks. The bifurcation diagrams of (a) the SOCS3-deleted network and (b) the SOCS3- and KLF4-deleted network, which fully agree with the prediction of SBA. Markers in different shapes denote distinct equilibrium branches as done in Fig. [5](#Fig5), where as the coloring for each chemical is determined by the biological interpretation (red for M1-like, and blue for M2-like). The vertical dashed lines denote the bifurcation threshold ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepa

Full text

Introduction

Various biological systems have been understood by employing mathematical approaches. Dynamical behaviors of biological systems are understood from interactions between biomolecules using mathematical models. For example, diversity of cell types and transition between them have been understood as bifurcation of dynamics of biological systems. Bifurcation is a phenomenon where the equilibria (or, more generally, the attractors) of a dynamical system undergo qualitative changes in response to continuous changes in system parameters.

Accumulated evidence has shown the connections between bifurcation theory and biological phenomena. For instance, Ozbudak et al. found that in Escherichia coli, the expression level of lactose permease (LacY), which is responsible for the cellular uptake of lactose and its analogues (e.g., TMG), switched significantly in response to continuous changes in the nutrient environment^1^. Ozbudak et al. developed a simple mathematical model and found that the jumps in LacY expression can arise from transitions between different steady states associated with the passage of a saddle-node bifurcation point. In this case, not only does the bifurcation create distinct cellular states, but it also creates a situation in which the transitions between these states are reversible, but the nutrient concentration thresholds that trigger the transitions are different in a state-dependent manner.

Another biological example associated with bifurcation is the polarization of macrophages. Macrophage is a type of white cell that can be polarized toward either an M1 (pro-inflammatory) or an M2 (anti-inflammatory) phenotype^2^. M1 and M2 macrophages can coexist under some microenvironmental conditions, which indicates that the dynamics of the polarization features the existence of alternative steady states^3^. Moreover, it has been known that a macrophage has the ability to switch its phenotype in response to stimuli^4–7^. These facts indicate that the dynamics of macrophage polarization allow transition between multiple steady states in response to parameter changes.

Bifurcation analyses of biological systems are mostly done in a heuristic manner. The analysis first starts from the construction of a mathematical model for the system under consideration, and then it is followed by the identification of equilibria in correspondence to given parameters. For a biological system, an analytical solution to determine the equilibrium point(s) by a function of parameters is not always technically available. Hence, to probe the bifurcation of a complex system, the equilibrium points are numerically determined while a parameter of the system is varying. This profiling step may capture bifurcation phenomena by chance; nonetheless, it is usually not known in advance that change in which of the numerous parameters can trigger a bifurcation. Also, it is not known in advance under what parameter condition bifurcation will occur, or which of the many variables will exhibit bifurcation behaviors.

The Jacobian matrix, denoted by J, plays a crucial role in the study of equilibrium bifurcations, as it takes the value of zero at a bifurcation point. Exploiting this property, a series of studies has provided fruitful findings that correlate the graphic representation of the network and the Jacobian determinant^8–10^. In addition, it has been shown that the number of alternative states can be even concluded from the sign of the Jacobian determinant^11,12^. While the Jacobian matrix can help identify the onset of a bifurcation, a numerical approach is still required unless the matrix is explicitly defined as a function of parameters. In practice, for large networked systems, conducting a bifurcation analysis remains challenging.

To cope with the practical problems of bifurcation analysis rooted in the application of the Jacobian determinant, Okada et al. developed Structural Bifurcation Analysis (SBA) by which the equilibrium bifurcation of a chemical reaction system is analyzed only from the structural information of the network^13,14^. Structural Bifurcation Analysis is based on a matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} determined by the network structure and a proved fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}=0$$\end{document} only when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J}=0$$\end{document} . This approach circumvents the problem that, without the functional form of reaction rate functions, the equilibria are not explicitly available for a bifurcation analysis. Nonetheless, a great deal about the possible bifurcation behaviors of the steady state can actually be determined from the network structure alone. Specifically, (1) By SBA, a network is decomposed into subnetworks based on a characteristics/index of sub-graphs, and this decomposition indicates all the conditions under which bifurcation can occur. (2) For each substructure given by the network decomposition, we can identify reactions containing the parameters that can induce the bifurcation. (3) Moreover, we can identify chemicals on the network that exhibit bifurcation behaviors when any substructure satisfies the bifurcation condition (Fig. 1).Fig. 1. Summary of Structural Bifurcation Analysis (modified from^13^). One can construct a matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} according to the network structure of a given chemical reaction networks. With the matrix, one can further conclude the conditions where bifurcations can occur, the inducing parameters, and the chemicals that exhibit the bifurcations.

That being said, there remained a couple of open problems in this approach: (i) the theorem was proved under the assumption that no conserved quantity is present in a given network, whereas the general case is open to discussion; (ii) the explicit relation between the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} and the Jacobian matrix is unclear, giving rise to a barrier to further development and applications. In this paper, we extend the theory of SBA by solving the two open problems. Our new theorem provides an explicit formula converting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} to a modified Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_g$$\end{document} (without redundant eigenvalues arising from the presence of conserved quantities, see Method). This breakthrough not only guarantees the application of SBA in general cases, but it also gives a corollary as a criterion for the (in-)stability of equilibrium points. Remarkably, our theorem is kinetics-independent and renders the study of system dynamics a hierarchical insight into network topology; furthermore, it provides a great potential to associate the framework of structural analysis with other theorems rooted in Jacobian matrices, e.g., multistationarity analysis.

Results

Unveiling bifurcation behaviors of chemical reaction networks through SBA

We begin by demonstrating how bifurcation behavior in chemical reaction networks can be predicted using SBA based on network topology. Suppose that a network composed of four reactions is given:

$$[eqn]$$

Here, the symbol \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\emptyset$$\end{document} denotes an external reservoir. The arrow “ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow$$\end{document} ” from chemical B to reaction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_2$$\end{document} in the network graph represents positive regulation. The subnetwork indicated by the red dashed block is a buffering structure identified by SBA, of which definition and properties are explained later. Then, according to the differences in the coefficients in the chemical reactions, we know the system follows the ordinary differential equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dfrac{d}{dt}\left( \begin{array}{c} x_A \\ x_B \end{array}\right) = \underbrace{\left( \begin{array}{cccc} 1 & -1 & 1 & -1 \\ 0 & 1 & -1 & 0 \end{array}\right) }_{\text {stoichiometric matrix } \nu } \left( \begin{array}{c} r_1(k_1) \\ r_2(x_A, x_B; k_2) \\ r_3(x_B; k_3) \\ r_4(x_A; k_4) \end{array}\right) {=}{:}f(\varvec{x}; \varvec{k}), \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_n$$\end{document} ’s are parameters of the reaction rate functions.

The SBA uses only the structural information from a reaction system, which is encapsulated in an augmented matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{A}}$$\end{document} (see Method Eq. (15) for the precise definition). In the absence of conserved quantities, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{A}}$$\end{document} is constructed based on the information whether reaction rates depend on each chemical (i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij} {:}{-} \frac{\partial{ r_i}}{\partial {x_j}} =0$$\end{document} or not) and a basis for the right kernel of the stoichiometric matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} . A kernel basis can be found from a perspective of graph theory by considering closed paths/loops, with the reservoir \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\emptyset$$\end{document} also taken as a node of the graph. In this case, we observe two loops: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_1: (\emptyset \xrightarrow [1]{\qquad} A \xrightarrow [2]{\qquad} B \xrightarrow [3]{\qquad} A \xrightarrow [4]{\qquad} \emptyset )$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_2: (A \xrightarrow [2]{\qquad} B \xrightarrow [3]{\qquad} A)$$\end{document} . Then, it is immediately followed by

$$[eqn]$$

By the theorem given by Okada, Tsai, and Mochizuki^13^, an equilibrium bifurcation occurs only when the determinant of the matrix takes the value of zero, since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}\propto \det {J_f}$$\end{document} in which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_f$$\end{document} is the Jacobian matrix of f in the governing equation Eq. (2). Therefore, the bifurcation points that can be identified by the condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J_f}=0$$\end{document} can also be identified by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {{\varvec{A}}}=0$$\end{document} .

One of the advantages of SBA over standard analysis using the Jacobian matrix is that it gives a network decomposition with the factorization of bifurcation conditions. A network can be decomposed when there exist buffering structures^15,16^, defined to be a subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} satisfying:

1. (Output-completeness) All the reactions regulated by any chemicals of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} are included in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} .
2. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 = \#(\text {chemicals in }\gamma )-\#(\text {reactions in }\gamma )+\#(\text {loops in }\gamma ) +\#(\text {conserved quantities in }\gamma )$$\end{document} . This is followed by a fact (which can be an alternative definition, see SI Appendix A) that the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} can be rearranged to be an upper block triangular matrix; that is, for some permutation matrices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_r$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_c$$\end{document} we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P_r \varvec{A} P_c = \left( \begin{array}{c|c} \varvec{A}_{\gamma } & *\\ \hline & \varvec{A}_{\gamma ^c} \end{array}\right) \end{aligned}$$\end{document}

such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}_\gamma$$\end{document} is a square block, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma ^c$$\end{document} denoting the complement of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} . The immediate consequence of the above factorization property is that a bifurcation occurs only if either (i) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_\gamma } = 0$$\end{document} or (ii) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\gamma ^c}} = 0$$\end{document} . For the considered example, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma {:}{=}\{B, R_2, R_3\}$$\end{document} is a buffering structure and one has

$$[eqn]$$

with the factorization of bifurcation conditions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \det {\varvec{A}} = \underbrace{\left| \begin{array}{cc} r_{2B} & 1 \\ r_{3B} & 1 \end{array}\right| }_{\det {\varvec{A}_\gamma }} \times \underbrace{\left| \begin{array}{cc} r_{4A} & 1 \\ 0 & 1 \end{array}\right| }_{\det {\varvec{A}_{\gamma ^c}}} = (r_{2B} - r_{3B})\cdot r_{4A}. \end{aligned}$$\end{document}

and it is clear that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {{\varvec{A}}_{\gamma ^{c}}}=r_{4A}>0$$\end{document} because A is a reactant of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_4$$\end{document} ; therefore, we can conclude that the system bifurcates only when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {{\varvec{A}}_{\gamma }}=r_{2B}-r_{3B}=0$$\end{document} .

Another advantage of SBA is its ability to predict bifurcation behavior without requiring kinetic information. With SBA^14^, it can be concluded that (i) when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_\gamma } = 0$$\end{document} , only chemicals in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} would exhibit the bifurcation behavior (SI Appendix A); (ii) if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\gamma ^c}} = 0$$\end{document} , all the chemicals could exhibit the bifurcation behavior. Nonetheless, in the considered example, we see that a bifurcation occurs only when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\gamma }} = r_{2B}-r_{3B} = 0$$\end{document} , and hence the behavior would only be exhibited by the chemical B. To verify that theoretical prediction holds regardless of kinetics, we plot the bifurcation diagrams of the common example network with three different reaction functions in ordinary differential equation models (SI Appendix B). As a parameter of the rate function of reaction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_2$$\end{document} is changing, the network exhibits three different types of bifurcations, whereas the bifurcation behaviors are always exhibited only by B as predicted (Fig. 2).Fig. 2. Bifurcation diagrams of the example network with different ODE models. The three models respectively show a (a) pitchfork, (b) transcritical, and (c) saddle-node bifurcation. Independent of the kinetics, we verify our prediction that the bifurcation behavior is exhibited only by B. The ODE models and the parameter settings are given by SI Appendix B.

The generalization of SBA with a standard procedure for chemical reaction networks

One of the challenges in bifurcation analyses in analytical approaches is posed when a system possesses conserved quantities. Conserved quantities create problems for standard bifurcation analysis because they cause the determinant of the standard Jacobian to vanish in such cases (as discussed in Method). Zero eigenvalues associated with these conserved quantities make \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J_f}\equiv 0$$\end{document} , rendering the determinant a poor indicator of bifurcation. Networks with conserved quantities have been studied in the structural analysis with an extended \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} with data on conserved quantities, but not in the context of bifurcation analyses but sensitivity analyses^16^. It was left as an open question whether such a structural approach still works for exploring bifurcation phenomena, that is, whether \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {{\varvec{A}}}$$\end{document} acts as a non-trivial bifurcation indicator despite the failure of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J_f}$$\end{document} .

Facing this problem, we find that the theory of SBA can be generalized to the case in which conserved quantities are present. Specifically, we prove that when conserved quantities are present, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} is proportional to a modified Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_g$$\end{document} that omits the redundant eigenvalues linked to conserved quantities. In other words, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} remains a reliable indicator of equilibrium bifurcations even in systems with conserved quantities (Method Eq. (18) and Eq. (19)). In particular, an explicit formula that transforms the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} into the modified Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_g$$\end{document} is given by our main theorem (Theorem 1 in Method), followed by an identity

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \det {\varvec{A}} = \dfrac{\det {(D^\intercal D)\cdot \det {(C^\intercal C)}}}{\det {\Lambda }}\cdot \det {J_g}, \end{aligned}$$\end{document}

where the constant matrices C, D, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} are structurally determined and defined in terms of the stoichiometric matrix. Furthermore, the system stability is assessed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {{\varvec{A}}}$$\end{document} with

Corollary 1

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\quad \mathop {\textrm{sign}}\nolimits (\det {\Lambda }\cdot \left. \det {\varvec{A}}\right| _{\varvec{x}}) \ne (-1)^{M-L} \implies \varvec{x}\text { is not a stable equilibrium point,}$$\end{document}

where M is the number of chemical species, and L is the number of conserved quantities (see Method).

In summary, our new mathematical findings demonstrate that the determinant of the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} serves as a general bifurcation and stability indicator for chemical reaction networks. This extends the applicability of SBA to a broader range of networks, including the problematic cases in which the standard Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_f$$\end{document} fails to be an informative indicator. Moreover, although the definition of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{A}}$$\end{document} differs in the presence of conserved quantities, the procedure suggested by reference^13^ for cases without conserved quantities can still be formally adopted without modification. To further streamline the workflow of SBA, we introduce a 6-step procedure (Fig. 3):Fig. 3. Summary of the standard procedure of Structural Bifurcation Analysis. The theory of SBA has been extended to networks with conserved quantities allowed, and we further improve the procedure for users to facilitate the roles and relations of subnetworks in the bifurcation dynamics.

1. Construct the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} according to the network topology.
2. Identify buffering structures and rearrange \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} into an upper block triangular form.
3. Decompose the network into subnetworks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} ’s and construct the influence graph according to the rearranged \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} .
4. For each subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} , compute \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\Gamma _i}}$$\end{document} .
5. Summarize the conditions in which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\Gamma _i}}$$\end{document} can change its sign for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} , and conclude the inducing parameters.
6. For each condition listed in (5), according to the influence graph, determine the bifurcating chemicals. We improve the protocol with the step of constructing the influence graph (defined in Method Eq. (28) and Eq. (29)), which provides an insight into the roles of and relations between subnetworks. With the generalization of the theory and the improvement of the procedure, in the following subsection, we would like to investigate the switching behaviors of a biological system at a cellular level.

Macrophage polarization network

Macrophages are specialized innate immune cells with significant plasticity to change their phenotypes in response to the microenvironment^17^. The phenotypic states of macrophages are broadly categorized into two distinct groups: M1 (pro-inflammatory) and M2 (anti-inflammatory), each wielding diverse functions in inflammation and tumor defense. M1 macrophages, characterized by highly activated STAT1 and NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B, produce pro-inflammatory cytokines and exert cytotoxic functions to combat pathogens^18,19^; on the other hand, M2, characterized by highly activated STAT3 and STAT6, represents an anti-inflammatory group which takes charge of tissue repairing and contributes to the resolution of inflammation^18,20,21^. In various diseases, the balance between M1 and M2 macrophages can be crucial and is thought to be strictly regulated according to the microenvironment including physiological and pathological conditions. Thus, for therapeutic purposes, researchers strive to comprehend and potentially manipulate macrophage polarization. For instance, it has been suggested that mutations in SOCS3 (suppressor of cytokine signaling 3), a gene that has a suppressive function for STAT3, may result in a hyperactivation of STAT3 and thus prevent the metastasis of skin cancer^22^. By elucidating the macrophage switching system, it would be possible to artificially switch macrophages between different phenotypes, which could potentially be used in the future for immunotherapy against cancer.

We take a simplified network associated with macrophage polarization as an example (Fig. 4). With the network, we first numerically verify that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} outperforms the standard Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_f$$\end{document} as a bifurcation indicator (by Theorem 1) and a stability indicator (by Corollary 1). After that, we demonstrate the power of SBA on predicting bifurcation behaviors based on the network topology. It is noteworthy that a signaling pathway as considered can usually possess conserved quantities, because signal transduction pathways often contain signaling proteins that reversibly transition between activated and inactivated states, and the total amounts of such (paired) proteins are constant in the timescale of chemical reactions. Specifically, in this macrophage polarization system, the total amount of each protein molecule (STAT1, NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B, STAT3, and STAT6) in a single cell corresponds to a quantity that does not change in a short period of time. Note the difference between the total amount and the activity level of each molecule. An important fact is that the expression levels of these proteins are known to vary among cells. In particular, it is reported that STAT1-deficiency in macrophage leads to a significantly low tendency of switching to an M1 phenotype^23,24^. Therefore, in the following analysis, we will focus on the qualitative changes in system behavior, or bifurcation behavior, that are caused by differences in expression levels between cells, i.e., differences in the total amount of signaling proteins (specifically, the total amount of inactivated and activated STAT1).Fig. 4. Macrophage polarization network for the wild type. By an arrow " \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dashv$$\end{document} " pointing from a complex to a reaction, we mean a negative regulation. The superscript \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+$$\end{document} denotes the phosphorylated/activated form for each of the four complexes STAT1, STAT3, STAT6, and NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B, of which signalings are essential to the macrophage M1/M2 polarization. In this simplified network, the M1 polarization is primarily associated with the dominance of STAT1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} and NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} , while the M2 polarization is associated with the dominance of STAT3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} and STAT6 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} .

SBA for the (wild-type) macrophage polarization network

The network of the macrophage polarization system, presented in Fig. 4 and referred to as the wild-type (WT) network, possesses no buffering structures (except for a trivial one, the entire network). Consequently, our structural method implies that the whole network is responsible for the bifurcation condition, and the determinant factor associated with the entire network, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} , is a bifurcation/stability indicator. Furthermore, it implies that if the system exhibits bifurcations, then all the four proteins must exhibit bifurcating behaviors.

To numerically verify these predictions from the SBA, we construct a model for the considered network. Reactions in the network can be categorized into two classes: deactivation (or degradation) and activation subject to negative regulation. For simplicity, the reaction rate function of deactivation is given by the law of mass action; on the other hand, the activation as well as the negative regulation involve phosphorylation, which is commonly modeled as a sigmoid function^25^, and this is considered in our model as well. The detailed parameter settings for the following numerical experiments can be found in supplementary materials (SI Appendix B).

For the convenience of the notation, the concentration of NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B is denoted by N, and that of activated form NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} is denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^p$$\end{document} . Similarly, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_i^p$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,3,6$$\end{document} ) denote the concentrations of inactivated and activated STAT family, respectively. Then, the four conserved quantities can be specified as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1 = S_1 + S_1^p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _2 = S_3 + S_3^p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _3 = S_6 + S_6^p$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _4 = N + N^p$$\end{document} . In order to examine the effect of variation in the amount of STAT1 expression, we numerically construct the bifurcation diagrams with varying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1$$\end{document} (Fig. 5 (a)). Specifically, for a fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1$$\end{document} , the equilibrium points are obtained by numerically finding the solution to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{0}=\nu \,\varvec{r}(\varvec{x})$$\end{document} . To numerically test the stability of each equilibrium point, say \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_0$$\end{document} , a small randomly given perturbation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \varvec{x}$$\end{document} such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\varvec{x}_0+\delta \varvec{x})$$\end{document} shares the same conserved quantities with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_0$$\end{document} is exerted to see whether the system state recovers to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_0$$\end{document} .

As shown in Fig. 5 (a), the numerical simulations demonstrate that all the chemicals exhibit the bifurcation behavior, which agrees with the SBA prediction. More specifically, there are two stable branches: one with highly-activated NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B and STAT1, and the other with highly-activated STAT3 and STAT6. Thus, our simple model successfully reproduced the known phenotypes of M1 and M2 macrophages. Besides, one of the stable branches, corresponding to M1 phenotype, disappears when STAT1 expression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1$$\end{document} decreases to a threshold value ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 6.58$$\end{document} ). This result is qualitatively consistent with the experimental observation that STAT1 deficiency in macrophages leads to a higher ratio of M2 to M1 polarization^23,24^.Fig. 5. Bifurcation analysis of the WT network. (a) The bifurcation diagrams with the conserved quantity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1=\text {STAT1}+\text {STAT1}^+$$\end{document} ranging in [6, 8.5] being the bifurcation parameter. Markers in different shapes are used to denote points on distinct equilibrium branches. The short and long stable equilibrium branches can be interpreted as the M1 (red) and M2 (blue) polarized states, respectively. (b) For each equilibrium as denoted in (a), we numerically estimate the determinants of the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} and the standard Jacobian. It is clear that the positiveness of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} represents the stability of equilibria, verifying our corollary of stability, and the standard Jacobian determinant always vanishes, failing to be an informative indicator of bifurcation. The vertical dashed line in each of panels illustrates the bifurcation threshold ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 6.58$$\end{document} ).

Collecting all the stable and unstable equilibrium points in the bifurcation diagram, we shall compute the corresponding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J_f}$$\end{document} to see whether the former does outperform \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J_f}$$\end{document} as an indicator for bifurcation and stability. Based on the network topology, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{A} = \left( {\begin{array}{cccccccc|cccc} r_{1S_1} & & & r_{1S_3^p} & & r_{1S_6^p} & & & 1 & & & \\ & r_{2S_1^p} & & & & & & & 1 & & & \\ & r_{3S_1^p} & r_{3S_3} & r_{3S_3^p} & & & & r_{3N^p} & & 1 & & \\ & & & r_{4S_3^p} & & & & & & 1 & & \\ & r_{5S_1^p} & & & r_{5S_6} & & & & & & 1 & \\ & & & & & r_{6S_6^p} & & & & & 1 & \\ & & & r_{7S_3^p} & & r_{7S_6^p} & r_{7N} & r_{7N^p} & & & & 1\\ & & & & & & & r_{8N^p} & & & & 1\\ \hline 1 & 1 & & & & & & & & & & \\ & & 1 & 1 & & & & & & & & \\ & & & & 1 & 1 & & & & & & \\ & & & & & & 1 & 1 & & & & \end{array}}\right) . \end{aligned}$$\end{document}

We numerically evaluate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {J_f}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} at all the equilibrium points and confirm that the system is unstable when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}\ge 0$$\end{document} , as claimed by Corollary 1 (Fig. 5 (b)). It is shown that the determinant of the standard Jacobian always vanishes and thus fails to be a decent bifurcation indicator; on the other hand, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}$$\end{document} stays informative of the occurrence of bifurcation.

Network alterations create buffering structures, resulting in localized bifurcation behaviors

As aforementioned, it has been suggested that the suppression of SOCS3 can cause the hyperactivation of STAT3 and prevent the metastasis of skin cancer^22^. Therefore, one may wish to examine the influence of deleting SOCS3 on the bifurcation behavior associated with macrophage polarization. SOCS3 can be induced by STAT1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} , STAT3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} , and NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} and negatively regulates STAT3 activation. In the WT network, presented in Fig. 4, these regulations are depicted by the three arrows acting on the activation of STAT3. In the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} corresponding to the WT network, these three regulations are expressed by negative partial derivatives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{3N^p}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{3S_1^p}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{3S_3^p}$$\end{document} . Therefore, the deletion of SOCS3 corresponds to removing these negative regulations from the network and setting these partial derivatives to zero in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} . As we will demonstrate below, this network alteration creates buffering structures, which can confine bifurcation behaviors.

We follow the 6-step standard procedure to predict the bifurcation behavior based on the network topology (see SI Appendix B for the altered network). In the first step, the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} is constructed as given in Eq. (8) but with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{3N^p} = r_{3S_1^p} = r_{3S_3^p} = 0$$\end{document} , since the regulation effects disappear due to the SOCS3 deletion. In step 2, the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} can be rearranged to be

$$[eqn]$$

According to the form, we define three subnetworks by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{array}{rcl} \Gamma _1 & =& \{\text {NF}\kappa \text {B}, \text {NF}\kappa \text {B}^+, R_7, R_8\} \\ \Gamma _2 & =& \{\text {STAT1}, \text {STAT1}^+, \text {STAT6}, \text {STAT6}^+, R_1, R_2, R_5, R_6\} \\ \Gamma _3 & =& \{\text {STAT3}, \text {STAT3}^+, R_3, R_4\}, \end{array} \end{aligned}$$\end{document}

Then, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} are two buffering structures, and the hierarchy of buffering structures with respect to the inclusion relations can be summarized as step 3 (Fig. 6). Next, calculating the determinants of square blocks in the diagonal in Eq. (9), then one obtains

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{array}{ccl} \det {\varvec{A}_{\Gamma _1}} & = & (r_{7N^p} - r_{7N} - r_{8N^p}) < 0 \\ \det {\varvec{A}_{\Gamma _2}} & = & (r_{1S_1} r_{5S_6} + r_{1S_1} r_{6S_6^p} - r_{1S_6^p} r_{5S_1^p} + r_{2S_1^p} r_{5S_6} + r_{2S_1^p} r_{6S_6^p}) \\ \det {\varvec{A}_{\Gamma _3}} & = & (r_{3S_3} + r_{4S_3^p})> 0. \end{array} \end{aligned}$$\end{document}

Notice that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{7N^p}<0$$\end{document} because NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} negatively regulates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_7$$\end{document} (i.e., the activation of NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B), and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{7N}>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{8N^p}>0$$\end{document} because the NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B and NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} are reactants in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_7$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_8$$\end{document} , respectively. This assures that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\Gamma _1}}$$\end{document} never becomes zero but stays negative, and we analogously have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\Gamma _3}}>0$$\end{document} . Therefore, as the conclusion of step 5, only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} is responsible for the occurrence of bifurcation. Finally, since the minimal buffering structure including \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} , the bifurcation behavior can be exhibited only in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} .Fig. 6. The influence graphs and the rearranged \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} matrices corresponding to the wild, the SOCS3-deleted, and the SOCS3-/KLF4- deleted types. In the rearranged matrix, we color the block \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\Gamma _2}$$\end{document} since the sign of its determinant can change, determining the occurrence of bifurcations. With the influence graphs, SBA predicts that chemicals in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} exhibit the bifurcation behavior in the SOCS3-deleted network, while the bifurcation behavior is limited within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} when KLF4-deletion is further applied to the network.

To verify the prediction of SBA that the SOCS3-deleted network can exhibit the bifurcation only in activated and inactivated forms of NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B, STAT1, and STAT6 (that is, chemicals in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} ), we adopt the ODE model built for the wild type but with the regulations mediated by SOCS3 removed (see SI Appendix B). The bifurcation diagram numerically computed aligns perfectly with the theoretical prediction based on SBA (Fig. 7 (a)). Furthermore, in our numerical analysis, the protein STAT3—located outside the buffering structure, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} , and therefore not exhibiting bifurcation—persists in a highly activated state, a characteristic akin to M2 polarization. This result is consistent with an experimental report that STAT3 is highly activated in the SOCS3-deleted mutants^22^.Fig. 7. Bifurcation analysis of the modified networks. The bifurcation diagrams of (a) the SOCS3-deleted network and (b) the SOCS3- and KLF4-deleted network, which fully agree with the prediction of SBA. Markers in different shapes denote distinct equilibrium branches as done in Fig. 5, where as the coloring for each chemical is determined by the biological interpretation (red for M1-like, and blue for M2-like). The vertical dashed lines denote the bifurcation threshold ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 7.25$$\end{document} ) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1$$\end{document} .

In the analysis of the SOCS3-deleted network described above, STAT3 is fixed in the M2-like phenotype due to a buffering structure that confines the bifurcation range within it. To further illustrate our method, we additionally modify this network to fix a subsystem to the M1-like phenotype. Equation (9) suggests that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} becomes an independent buffering structure if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{7S^p_6}\equiv 0$$\end{document} in the SOCS3-deleted network. This partial derivative is corresponding to the regulation of STAT6 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} on the activation of NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B via the gene expression of KLF4 (Krüppel-like factor 4). When the deletions of both SOCS3 and KLF4 are applied (see SI Appendix B for the doubly-altered network), we still have Eq. (11); namely, only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} is responsible for the occurrence of bifurcation. Moreover, since the minimal buffering structure including \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} itself, SBA concludes that only the four chemicals STAT1, STAT1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} , STAT6, and STAT6 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^+$$\end{document} can exhibit the bifurcation behaviors. With an analogous numerical simulation, we once again verify that the prediction of SBA (Fig. 7 (b)). Only STAT1 and STAT6, located within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} , exhibit bifurcation, resulting in branches where both of the proteins are M1-like or both are M2-like. Conversely, STAT3 and NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B, located outside \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} , remain persistently highly activated, each resembling the M2 and M1 phenotypes, respectively. In other words, our results suggest that the double mutations would create a phenotype that simultaneously displays M1- and M2-associated features, robustly (e.g. independently of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1$$\end{document} ).

Discussion

In this paper, we further extend structural bifurcation analysis, which determines the fundamental properties of bifurcation behavior in chemical reaction systems based solely on network structures. We have increased the theoretical generality of the method while also demonstrating its practicality for elucidating biological systems. We have shown, through rigorous mathematical formulae, the direct relationship between the Jacobian matrix, which characterizes conventional theory, and the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} , which characterizes structural bifurcation analysis. This not only clarifies the relationship with conventional theory, but also makes it possible to analyze bifurcation behavior in response to changes in conserved quantities in addition to changes in reaction parameters. We also showed the method of explicitly obtaining the influence relationship of parameter perturbation on substructures on the network as an influence graph from the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} . Furthermore, we showed the algorithm for comprehensively performing structural bifurcation analysis and structural sensitivity analysis. These results were applied to the macrophage polarization system, and not only did they clarify the mechanisms behind cell fate regulation, but they also enabled us to predict changes in cell fate switching behavior in mutants.

As the example of the macrophage polarization system shows, bifurcation has important significance for biological systems, as it creates discrete states and transitions between them. Nevertheless, bifurcation analysis has not yet achieved sufficient popularity in the life sciences. The main reason for this is the difficulty of the analysis. To find bifurcation phenomena, it is necessary to identify the equilibrium points of the target system and investigate the dependence of the equilibrium points on parameters. To perform these analyses on biological systems, which include many molecular species, i.e., high-dimensional dynamical systems, the practical way is to perform numerical analysis. Nonetheless, such a high-dimensional biological system always comes with a large number of parameters involved. Therefore, until now, the only way to identify the equilibrium points and parameter dependencies of high-dimensional dynamical systems in life was through a trial-and-error approach using numerical calculations.

The structural bifurcation analysis we propose has the potential to address the aforementioned issues by determining fundamental bifurcation properties solely based on network structure in a model-free manner. This approach becomes especially powerful when buffering structures, identified according to the network structures, are present. By decomposing large biological systems into substructures determined by buffering structures, we can elucidate the structural conditions necessary for bifurcation, and pinpoint the bifurcation parameters and variables that exhibit bifurcation behavior within each substructure of the network.

The matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} directly reflects the topology of the network and provides deeper insight into the network structure for analysis. The new mathematical finding presented in this paper not only generalizes the theory of structural bifurcation analysis to chemical reaction networks with conserved quantities but also opens avenues for further expansion of the theory. For instance, it implies that the (in-)stability of the system can be inferred from the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} . We anticipate that this breakthrough can be integrated into other methods rooted in the analyses of the Jacobian matrix. For example, the type of bifurcation can be determined by the eigenspectrum of the Jacobian matrix^26^, and the sign of the Jacobian determinant can give a hint whether the system possesses a single equilibrium point^27^. We believe that using the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} instead of the Jacobian matrix can render these theories a structural insight into chemical reaction networks.

With the network in Example 1, we obtained three different bifurcation scenarios by assuming three different mathematical models, i.e., sets of reaction rate functions. While different reaction rate functions yielded different bifurcation types, the use of a common network produced consistent results for bifurcation conditions, bifurcation parameters, and bifurcating variables. This clearly demonstrates that these three dynamical properties are determined solely by the structure of the reaction network. Historically, bifurcation patterns have been a significant focus in the mathematical sciences. However, within the context of structural bifurcation analysis, bifurcation patterns may be regarded as incidental properties that emerge only after the network structure has been defined and a specific model is assumed.

In this paper, we focus on equilibrium bifurcations, which are characterized by the Jacobian matrix acquiring a zero eigenvalue. Our analysis does not extend to bifurcations involving complex eigenvalues, such as Hopf bifurcations, whose onset is indicated by the emergence of purely imaginary eigenvalues. Such eigenvalue configurations are not captured by the determinant of the Jacobian, and the framework of SBA is thus not applicable in these cases. As shown in Okada et al.^14^, in the absence of conserved quantities, the characteristic polynomial of the Jacobian matrix J can be written in terms of the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} . Our main theoretical result (Theorem 1 in Methods), which connects the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} to the modified Jacobian—a system without conserved quantities—therefore implies that the full eigenvalues can be examined structurally even in the presence of conserved quantities. We consider the study of bifurcations associated with complex eigenvalues to be an important and promising avenue for future research.

The structural bifurcation analysis presented in this study will serve as a powerful tool for advancing our understanding of biological systems. To demonstrate this, we analyzed the macrophage polarization system and derived general predictions from the network structures. In wild-type cells, macrophage polarization was examined by varying the expression level of STAT1. When STAT1 expression is low, macrophages are predicted to exhibit only the M2 type (anti-inflammatory), lacking diversity in phenotype. Conversely, when STAT1 expression is high, both M1 (inflammatory) and M2 types appear. In SOCS3 mutants, our analysis indicates that STAT3 keeps highly activated, independent of STAT1 expression, without exhibiting bifurcation behavior. This finding aligns with experimental observations of STAT3 hyperactivation in SOCS3 mutants. However, it is important to note that bifurcation behavior can still occur in the activities of STAT1 and STAT6 in SOCS3 mutants. This result suggests that the polarization behavior of macrophages with SOCS3 deletion may differ depending on which gene’s expression is observed. For instance, STAT3 may always display hyperactivation and an M2 phenotype, while STAT1 and STAT6 show normal polarization into either M1 or M2, similar to the wild type. Additionally, structural bifurcation analysis suggests that the suppression or deletion of KLF4 in the SOCS3-knockout system further eliminates multistability in both NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B and STAT3, leading to their simultaneous hyperactivation. Interestingly, while NF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa$$\end{document} B hyperactivation is associated with the induction of pro-inflammatory genes, STAT3 hyperactivation promotes anti-inflammatory responses. This suggests that a macrophage lacking both SOCS3 and KLF4 may exhibit seemingly contradictory behaviors, warranting further investigation to understand its role in inflammation.

Actually, the reliability of network information in databases is widely concerned. Networks may contain unknown reactions or regulatory interactions that are not reflected in existing databases. Our structural approach is particularly powerful in such challenging circumstances, precisely because SBA relies solely on the network topology. Suppose that SBA is applied to a network on a database, but experimental results in the verification phase do not align with the prediction. In that case, misinformation about the network structure presented in the databases can be firmly concluded. In the subsequent phase of network refinement, one can search for missing reactions or regulations whose inclusion would reconcile the predictions of SBA with experimental observations. By experimentally testing these candidates, the true underlying network structure can be elucidated.

We believe that our predictions on the bifurcation behaviors of macrophage polarization can be verified through experiments and that these experiments can provide valuable results. As noted, experimental validations can help identify unknown reactions/regulations and guide network refinement. Finally, it is noteworthy that the M1 and M2 dichotomy is a simplistic and convenient categorization of possible macrophage phenotypes. There are subtypes between the M1 and M2 phenotypes, and their immunological roles are increasingly recognized as important^28,29^. Our structural method may also be beneficial in predicting potential subtypes based on a network structure or in proposing modifications to achieve targeted phenotypes.

Methods

Review of Structural Bifurcation Analysis (SBA)

Suppose that a network is given by N reactions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} R_n: \quad \nu ^-_{1,n}X_1+\ldots +\nu ^-_{M,n}X_M \xrightarrow {\hspace{1cm}} \nu ^+_{1,n}X_1+\ldots +\nu ^+_{M,n}X_M, \qquad (n=1,\ldots ,N) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{X_m\}_{m=1}^M$$\end{document} is the collection of chemical species and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu ^\pm _{m,n}$$\end{document} ’s represent the coefficients in the reactions. Then, an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M\times N$$\end{document} stoichiometric matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} is defined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\nu _{m,n}) {:}{=}(\nu ^+_{m,n} - \nu ^-_{m,n})$$\end{document} , encoding the network topology as a weighted incidence matrix of the directed graph. Then, the dynamics follows the ordinary differential equations

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \displaystyle \frac{d}{dt}\,x_m = \sum _{n=1}^{N} \nu _{m,n} r_{n}(\varvec{x};\varvec{k}) {=}{:}f_m(\varvec{x};\varvec{k}) \end{aligned}$$\end{document}

for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m=1, 2,\ldots ,M$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{k}$$\end{document} is a parameter vector in the reaction rate function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{r}:\mathbb {R}^M\times \mathbb {R}^N\rightarrow \mathbb {R}^N$$\end{document} . There is no universal rule of determining the functional form for reaction rate function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{r}$$\end{document} . The nonlinearity of reaction rate functions in general contributes to the difficulties of modeling approaches in quantitatively estimating the parameters based on data. The influences of chemicals to reaction rates, meanwhile, can usually be qualitatively determined by the biological background knowledge. For instance, the value of the partial derivative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial {r_n}}{\partial {x_m}}$$\end{document} is positive if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_m$$\end{document} is a reactant in the reaction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_n$$\end{document} , and it is negative if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_m$$\end{document} negatively regulates the reaction as a suppressor, and it can also be zero when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_m$$\end{document} does not affect the reaction rate of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_n$$\end{document} .

Without functional forms of kinetics, the equilibrium dynamics can still be well characterized to a degree with the help of two linear subspaces. When the system is at an equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\varvec{x}}\in \mathbb {R}^M$$\end{document} , the change rate is zero and thus we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{r}(\overline{\varvec{x}}) = \sum _{p=1}^{P} \mu _p \varvec{c}^p \end{aligned}$$\end{document}

for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \in \mathbb {R}^P$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\varvec{c}^p\}_{p=1}^{P}$$\end{document} is a basis of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ker {\nu }\subseteq \mathbb {R}^N$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P=\dim \ker {\nu }$$\end{document} . We may put \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C{:}{=}(\varvec{c}^1, \ldots , \varvec{c}^P)\in \mathbb {R}^{N\times P}$$\end{document} and let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle C\right\rangle$$\end{document} denote the space spanned by column vectors of C. Besides, we say a chemical reaction network has a conserved quantity if a linear combination of chemicals is time-invariant, namely, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\equiv \varvec{d}^\intercal \dot{\varvec{x}}$$\end{document} for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{d}\in \mathbb {R}^M$$\end{document} for any \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}\in \mathbb {R}^M$$\end{document} . By the dynamics given by Eq. (13), we see \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{d}\in \ker {\nu ^\intercal }\simeq \mathop {\textrm{coker}}\nolimits {\nu }$$\end{document} . Therefore, take a cokernel basis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D = (\varvec{d}^1,\ldots ,\varvec{d}^Q)\in \mathbb {R}^{M\times Q}$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q=\dim \ker {\nu ^\intercal }$$\end{document} , then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\eta }{:}{=}D^\intercal \varvec{x}$$\end{document} specifies all the conserved quantities. We similarly denote \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle D\right\rangle =\ker {\nu ^\intercal }$$\end{document} .

With the graph and the biological information rendered by the chemical reactions, a square matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} is defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{A} {:}{-} \left( \begin{array}{c|c} \begin{array}{ccc} {} & {} & {} \\ {} & \dfrac{\partial \varvec{r}}{\partial \varvec{x}} & {} \\ {} & {} & {} \end{array} & C \\ \hline D^\intercal & \end{array}\right) , \end{aligned}$$\end{document}

of which lower-right block is empty (throughout the article, trivial blocks are denoted empty). Together with C indicating loops of reactions, D determining the conserved quantities, the zero-nonzero pattern of entries of the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} directly encodes the structure of the reaction network structure. Such a matrix was first introduced for the sensitivity analysis of chemical reaction networks^30,31^; specifically, a sensitivity matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S} = -\varvec{A}^{-1}$$\end{document} determines the system responses to perturbations in parameters. More explicitly,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \dfrac{\delta {x_i}}{\delta {k_j}} = (\varvec{S})_{ij} \end{aligned}$$\end{document}

indicates the change \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x_i\mapsto x_i+\delta x_i)$$\end{document} in the chemical \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_i$$\end{document} in response to a parameter change \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k_i\mapsto k_i+\delta k_i)$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j$$\end{document} (see SI Appendix A.I for details of structural sensitivity analysis).

The identification of buffering structures is important in the approach of structural analysis, since a localization principle arises from buffering structures—the influence of a parameter change inside a buffering structure localizes inside the buffering structure (see SI Appendix A.III.). The localization principle does not hold only in the context of the aforementioned sensitivity problem, but it was extended to bifurcation analysis by the mathematical finding of Okada et al.^13,14^ that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \varvec{A}\propto \det {J_f}$$\end{document} when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} is full-ranked (i.e., the system has no conserved quantities), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_f$$\end{document} being the Jacobian of f in the Eq. (13). This leads to the development of Structural Bifurcation Analysis, in which a bifurcation is identified by the vanishing of the determinant of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} . In particular, as shown in Eq. (4), the identification of buffering structures facilitates the factorization of the determinant (followed by the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \varvec{A}=0\iff$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \varvec{A}_\gamma =0$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \varvec{A}_{\gamma ^c}=0$$\end{document} ). The presence of multiple buffering structures (say \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _\alpha$$\end{document} ’s) highly benefits us in bifurcation analysis with the capability of

1. Factorizing the bifurcation conditions: the bifurcation analysis of a complex network can be facilitated by the network decomposition into determinant structures (denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} ’s), which are defined as buffering structures with subtraction of their inner buffering structures (Fig. 1 (3)). For each determinant structure, the condition for bifurcation occurrence can be examined from its network structure.
2. Finding (bifurcation-)inducing parameters: for each determinant structure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} , the bifurcation parameters associated with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \varvec{A}_{\Gamma _{\beta }}=0$$\end{document} , termed as the inducing parameters for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} , are identifiable on the network. The bifurcations are induced by the parameters which are outside any buffering structures non-intersecting with the determinant structure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} (Fig. 1 (4)).
3. Determining the bifurcating chemicals: for each determinant structure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} , the chemicals exhibiting bifurcation behaviors associated with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det \varvec{A}_{\Gamma _{\beta }}=0$$\end{document} , termed as the bifurcating chemicals for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} , are identifiable on the network. When the bifurcation condition is satisfied, the bifurcation of steady-state concentrations (and fluxes) appears only inside the (minimal) buffering structure that contains the determinant structure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _{\beta }$$\end{document} (Fig. 1 (5)).

The generalization of SBA

SBA is rooted in the finding of Okada et al. that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}}\propto \det {J_f}$$\end{document} , and we remark that the theorem is proved under the assumption that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} is full-ranked, that is, no conserved quantity is present in the considered network^13,14^. It remains an open problem whether the framework of SBA can be employed when conserved quantities exist.

What is the problem with the conserved quantities? When the cokernel space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle D\right\rangle$$\end{document} is nontrivial, we have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D^\intercal J_f = D^\intercal \nu \frac{\partial \varvec{r}}{\partial \varvec{x}} = (\nu ^\intercal D)^\intercal \frac{\partial \varvec{r}}{\partial \varvec{x}} \equiv \varvec{0}$$\end{document} . This implies that the Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_f$$\end{document} is always non-invertible because of the vanishing eigenvalues due to the nontrivial cokernel space. Intuitively speaking, a trajectory \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}(t)$$\end{document} always lies on a hyperplane parallel to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {\textrm{im}}\nolimits \nu$$\end{document} , which does not span the whole space when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} is not full-ranked. Therefore, L of the eigenvalues that measure diverging rates of trajectories must be zero, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L {:}{-} \dim \mathop {\textrm{coker}}\nolimits \nu = (M-\dim \mathop {\textrm{im}}\nolimits \nu )$$\end{document} is the number of independent conserved quantities.

To circumvent the problem with the redundant eigenvalues, the intuition above motivates us to have a change of coordinates with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\mathop {\textrm{im}}\nolimits \nu )$$\end{document} and its orthogonal complement \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\mathop {\textrm{im}}\nolimits \nu )^\perp$$\end{document} being the two axes. To do so, we put V to be a full-ranked matrix such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle V\right\rangle =\mathop {\textrm{im}}\nolimits \nu$$\end{document} and consider the pseudo-inverses

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} V^\dagger {:}{=}(V^\intercal V)^{-1} V^\intercal \quad \text {and}\quad D^\dagger {:}{=}(D^\intercal D)^{-1} D^\intercal . \end{aligned}$$\end{document}

Then, a linear isomorphism for the desired change of coordinate is obtained as

$$[eqn]$$

or, alternatively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}=V\varvec{\xi } + (D^\dagger )^\intercal \varvec{\eta }$$\end{document} . Then, the dynamics can be rewritten as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{array}{cccc} \dot{\varvec{\xi }} & = & g(\varvec{\xi };\varvec{\eta },\varvec{k}) & {:}{-} V^\dagger \nu \,\varvec{r}(V\varvec{\xi } + (D^\dagger )^\intercal \varvec{\eta };\varvec{k}) \\ \dot{\varvec{\eta }} & = & \varvec{0} \end{array}\right. . \end{aligned}$$\end{document}

This gives a clear picture showing that the dynamics of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}$$\end{document} solely counts on that of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }$$\end{document} , and thus it is sufficient to consider the Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_g=V^\dagger {J_f}V$$\end{document} of g for the bifurcation and stability analyses. By considering \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_g$$\end{document} instead of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_f$$\end{document} , the problem of the vanishing eigenvalues due to conserved quantities is well circumvented. Moreover, the Jacobian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_g$$\end{document} is associated with the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} which is given by the network topology. Considering a constant matrix

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Lambda {:}{-} \left( \begin{array}{c|c} \begin{array}{ccc} {} & {} & {} \\ {} & \nu & {} \\ {} & {} & {} \end{array} & D \\ \hline C^\intercal & \end{array}\right) , \end{aligned}$$\end{document}

we have the relation associating \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} and the Jacobian, given by

Theorem 1

Given a kernel basis C and a cokernel basis D of the mapping \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} , then there exists an orthogonal basis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{T}\in \mathbb {R}^{(M+Q)\times (M+Q)}$$\end{document} such that and hence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}} = \dfrac{\det {(D^\intercal D)\cdot \det {(C^\intercal C)}}}{\det {\Lambda }}\cdot \det {J_g}$$\end{document} .

In particular, when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dim \mathop {\textrm{coker}}\nolimits \nu =0$$\end{document} , we have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}} = \dfrac{\det {(C^\intercal C)}}{\det {\Lambda }}\cdot \det {J_f}$$\end{document} .

The relation further gives Corollary 1 as a criterion for stable equilibria. The proofs are left in SI Appendix C.

An example network with conserved quantities

To show the problem induced by the presence of conserved quantities and how our theorem circumvent it, we consider a simple network.

$$[eqn]$$

The corresponding stoichiometric matrix is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nu = \left( {\begin{array}{cccccc} 1 & -1 & {\,} & {\,} & {\,} & {\,} \\ {\,} & 1 & -1 & {\,} & 1 & -1 \\ {\,} & {\,} & 1 & -1 & {\,} & {\,} \\ {\,} & {\,} & {\,} & {\,} & -1 & 1 \\ {\,} & {\,} & {\,} & {\,} & -1 & 1 \\ \end{array}}\right) . \end{aligned}$$\end{document}

Furthermore, one sees that the cokernel space is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ker \nu ^\intercal = \langle \,(0,0,0,1,-1)^\intercal \,\rangle , \end{aligned}$$\end{document}

since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _1 {:}{-} 1\cdot {x_D} + (-1)\cdot {x_E}$$\end{document} determines the conserved quantity. On the other hand, the kernel basis can be obtained by considering the two loops \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_1: (B \xrightarrow [6]{\qquad} (D+E) \xrightarrow [5]{\qquad} B)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_2: (\emptyset \xrightarrow [1]{\qquad} A \xrightarrow [2]{\qquad} B \xrightarrow [3]{\qquad} C \xrightarrow [4]{\qquad} \emptyset )$$\end{document} . As such, one obtains

$$[eqn]$$

By the identification as given in SI Appendix A, it can be seen that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\Gamma _1\cup \Gamma _2)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _3$$\end{document} are three buffering structures, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1 {:}{-} \{A, R_2\}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2{:}{-}\{C, R_4\}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _3 {:}{-} \{B, D, E, R_3, R_5, R_6\}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _4 {:}{-} \{R_1\}$$\end{document} . We have the rearranged form

$$[eqn]$$

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_r$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_c$$\end{document} are permutation matrices with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {P_r} = \det {P_c} = -1$$\end{document} . The step of rearranging \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} according to the buffering structure highly facilitates the computation and factorization of the determinants, giving

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \det {\varvec{A}} = \det {P_r\varvec{A}P_c} = \prod _{i=0}^3 \det {\varvec{A}_{\Gamma _i}} = r_{2A}\,r_{4C}\,(r_{3B}r_{5D}+r_{3B}r_{5E}+r_{3E}r_{6B}), \end{aligned}$$\end{document}

in which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma_0 {:}{-} \left( \bigcup _{i=1}^3\Gamma _i\right) ^c$$\end{document} is the complement of the union of all buffering structures. With Theorem 1, it is known that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}} = \frac{\det {(D^\intercal D)\cdot \det {(C^\intercal C)}}}{\det {\Lambda }}\cdot \det {J_g} = \det {J_g}$$\end{document} since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {(C^\intercal C)} = 8$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {(D^\intercal D)} = 2$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\Lambda } = 16$$\end{document} ; meanwhile,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \det {J_f} = \left| { \begin{array}{ccccc} -r_{2A} & & -r_{2C} & & \\ r_{2A} & -r_{3B} - r_{6B} & r_{2C} & r_{5D} & -r_{3E} + r_{5E} \\ & r_{3B} & -r_{4C} & & r_{3E} \\ & r_{6B} & & -r_{5D} & -r_{5E} \\ & r_{6B} & & -r_{5D} & -r_{5E} \end{array}}\right| \equiv 0, \end{aligned}$$\end{document}

which always vanishes and thus fails to be an informative bifurcation indicator.

Despite the vanishing standard Jacobian determinant, we see that the generalized SBA can be applied to this network to predict the bifurcation behaviors. Since A is a reactant in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_2$$\end{document} and C is a reactant in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_4$$\end{document} , both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{2A}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{4C}$$\end{document} do not vanish but stay positive. The Eq. (26) implies that a bifurcation occurs only when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det {\varvec{A}_{\Gamma _3}}=0$$\end{document} . Furthermore, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _3$$\end{document} is a buffering structure, one can conclude that the bifurcation behaviors can be exhibited only by chemicals B, D, and E.

Influence graph among subnetworks

As shown by the two pedagogical examples, a network decomposition into subnetworks readily follows the rearrangement of the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{A}$$\end{document} . Explicitly, the square blocks in the diagonal of the upper block triangular \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_r\varvec{A}P_c$$\end{document} define disjoint subnetworks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} ’s. More specifically, the rearranged matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_r\varvec{A}P_c$$\end{document} takes the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P_r \varvec{A} P_c = \left( \begin{array}{cccc} A_{\Gamma _1} & Y_{12} & \cdots & Y_{1S}\\ {\,} & A_{\Gamma _2} & \cdots & Y_{2S} \\ {\,} & {\,} & \ddots & \vdots \\ {\,} & {\,} & {\,} & A_{\Gamma _S} \end{array}\right) , \end{aligned}$$\end{document}

and the subnetworks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} ’s can be determined accordingly: a subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} consists of a chemical/reaction if a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} ’s column/row is corresponding to the chemical/reaction. We term the subnetworks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _\alpha$$\end{document} ’s the determinant structures, which are not necessarily buffering structures, but a buffering structure is always a union of some determinant structures. In other words, determinant structures give a finer network decomposition then buffering structures do.

Our interest is the influence relation among these subnetworks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} ’s, which is termed determinant structures. Before introducing how the relation is revealed by SBA, we shall first clarify what it means by influence. We say a subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _\alpha$$\end{document} influence another subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _\beta$$\end{document} if there exists some parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_i\in \Gamma _\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_j\in \Gamma _\beta$$\end{document} such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\delta x_i}{\delta k_j}\ne 0$$\end{document} , of which idea is identical to Eq. (16).

The influence relation among the determinant structures is encoded by blocks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{ij}$$\end{document} ’s in the matrix representation (28). In particular, a subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _j$$\end{document} influences another subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{ij}$$\end{document} is not a zero matrix. As such, we can draw a directed graph to visualize the influencing relation among subnetworks, by considering the adjacency matrix G^32^ defined by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} G_{ij} = \left\{ \begin{array}{cl} 1 & \text {if }i\ne j\text { and }Y_{ji}\text { is not a zero block}, \\ 0 & \text {otherwise.} \end{array}\right. \end{aligned}$$\end{document}

We may take the network given in Eq. (21) for example, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1 {:}{=}\{A, R_2\}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2{:}{=}\{C, R_4\}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _3 {:}{=}\{B, D, E, R_3, R_5, R_6\}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _4{:}{=}\{R_1\}$$\end{document} , according to the rearranged form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_r\varvec{A}P_c$$\end{document} as in Eq. (25). The matrix can be rewritten into the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P_r\varvec{A}P_c = \left( \begin{array}{cccc} \varvec{A}_{\Gamma _1} & Y_{12} & {\,} & Y_{14} \\ {\,} & \varvec{A}_{\Gamma _2} & {\,} & Y_{24} \\ {\,} & {\,} & \varvec{A}_{\Gamma _3} & Y_{34} \\ {\,} & {\,} & {\,} & \varvec{A}_{\Gamma _4} \end{array}\right) , \end{aligned}$$\end{document}

suggesting the adjacency matrix G of the influence graph as

$$[eqn]$$

so that the influence relation among determinant structures can be clearly visualized. A determinant structure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} influences another subnetwork \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _j$$\end{document} if there is a directed path from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _i$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _j$$\end{document} in the influence graph. With the graph, it is straightforward to see that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1, (\Gamma _1\cup \Gamma _2)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _3$$\end{document} are buffering structures, since they do not have their downstreams in the influence graph. It is clear to see that parameter changes in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} could possibly influence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1$$\end{document} but not \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _3$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _4$$\end{document} , since that only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _1$$\end{document} is in the downstream for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma _2$$\end{document} . It is also worthy of mention that the graph is acyclic (i.e., with no directed cycles), which means that a subnetwork would not be simultaneously in the upstream and the downstream for any other subnetwork. In summary, the influence graph renders the network a hierarchy of subnetworks, giving a structural insight into the equilibrium dynamics.

Supplementary Information

Supplementary Information.