Sensitivity and Network Topology in Chemical Reaction Systems

TL;DR

This paper extends a framework for analyzing biochemical reaction networks to include conserved concentrations, introduces new measures of robustness, and demonstrates that E. coli's metabolic network is more robust than random networks.

Contribution

It generalizes the law of localization to reaction systems with conserved concentrations and proposes novel network metrics for robustness analysis.

Findings

E. coli metabolic network exhibits higher robustness than random networks.

The generalized framework applies to any reaction system with conserved quantities.

Network topology influences response patterns and robustness.

Abstract

In living cells, biochemical reactions are catalyzed by specific enzymes and connect to one another by sharing substrates and products, forming complex networks. In our previous studies, we established a framework determining the responses to enzyme perturbations only from network topology, and then proved a theorem, called the law of localization, explaining response patterns in terms of network topology. In this paper, we generalize these results to reaction networks with conserved concentrations, which allows us to study any reaction systems. We also propose novel network characteristics quantifying robustness. We compare E. coli metabolic network with randomly rewired networks, and find that the robustness of the E. coli network is significantly higher than that of the random networks.

Equations

$$\frac{d{x}_m}{dt}=\sum _{j=1}^R{\nu }_{mj}{r}_j(k_j;x).$$

$$\displaystyle j:\ \sum_{m=1}^{M}y^{j}_{m}X_{m}\rightarrow\sum_{m=1}^{M}\bar{y}^{j}_{m}X_{m}\ \$$

$$(X_m:\mathrm{types}\mathrm{of}\mathrm{molecules}),$$

$${\nu }_{mj}\equiv {\bar{y}}_m^j-y_m^j.$$

$$\frac{\partial r_j}{\partial x_m}>0\mathrm{if}{y}_m^j>0,\mathrm{otherwise}\frac{\partial r_j}{\partial x_m}=0.$$

$$r_j=\sum _{\alpha =1}^N{\mu }^{\alpha }{c}_j^{\alpha }$$

$$\displaystyle l_{a}\equiv{\mbox{\boldmath$d$}_{a}}\cdot{\mbox{\boldmath$x$}}$$

$$\{K_1,\dots ,K_R,K_{R+1},\dots ,K_{R+N_c}\}$$

$$\equiv \{k_1,\dots ,k_R,l_1,\dots ,l_{N_c}\}.$$

$$\displaystyle{\bf A}\left(\begin{array}[]{cccc}\delta_{1}{\mbox{\boldmath$x$}}&\ldots&\delta_{R+N_{c}}{\mbox{\boldmath$x$}}\\ \hline\cr\delta_{1}{\mbox{\boldmath$\mu$}}&\ldots&\delta_{R+N_{c}}{\mbox{\boldmath$\mu$}}\\ \end{array}\right)=-\left(\begin{array}[]{c| c}{\bf E}_{R}&{\bf 0}\\ \hline\cr{\bf 0}&{\bf E^{\prime}}_{N_{c}}\end{array}\right),$$

$$\displaystyle\delta_{J^{*}}\mbox{\boldmath$x$}=(\delta_{J^{*}}x_{1},\ldots,\delta_{J^{*}}x_{M})^{t},$$

$$\displaystyle\delta_{J^{*}}\mbox{\boldmath$\mu$}=(\delta_{J^{*}}\mu^{1},\ldots,\delta_{J^{*}}\mu^{N})^{t},$$

$$\displaystyle{\bf A}\equiv\begin{array}[]{c}\\ \ _{R}\Bigg{\updownarrow}\\ \\ \ _{N_{c}}\Bigg{\updownarrow}\\ \\ \end{array}\underset{\underset{M}{\ext@arrow 0099\arrowfill@\longleftarrow\relbar\longrightarrow{}{\hskip 40.00006pt}}\ \ \ \ \ \ \underset{N}{\ext@arrow 0099\arrowfill@\longleftarrow\relbar\longrightarrow{}{\hskip 50.00008pt}}}{\left(\begin{array}[]{cccc|c }&&&&\\ &\mbox{\smash{ $r_{jm}$}}&&&{-{\mbox{\boldmath$c$}}^{\,1}\ {\ldots}\ -{\mbox{\boldmath$c$}}^{\,N}}\\ &&&&\\ \hline\cr&{-({\mbox{\boldmath$d$}}^{\,1})^{T}}&&&\\ &\vdots&&&\mbox{\smash{\ $\bf 0$}}\\ &{-({\mbox{\boldmath$d$}}^{\,N_{c}})^{T}}&&&\end{array}\right)},$$

$$\displaystyle\left(\begin{array}[]{cccc}\delta_{1}{\mbox{\boldmath$x$}}&\ldots&\delta_{R+N_{c}}{\mbox{\boldmath$x$}}\\ \hline\cr\delta_{1}{\mbox{\boldmath$\mu$}}&\ldots&\delta_{R+N_{c}}{\mbox{\boldmath$\mu$}}\\ \end{array}\right)=-{\bf A}^{-1}\left(\begin{array}[]{c| c}{\bf E}_{R}&{\bf 0}\\ \hline\cr{\bf 0}&{\bf E^{\prime}}_{N_{c}}\end{array}\right).$$

$${\delta }_{J^{\ast }}{r}_j=\sum _{\alpha =1}^N{\delta }_{J^{\ast }}{\mu }^{\alpha }{c}_j^{\alpha },$$

$$\displaystyle\left(\begin{array}[]{c c c}\delta_{1}{\mbox{\boldmath$r$}}&\ldots&\delta_{R+N_{c}}{\mbox{\boldmath$r$}}\end{array}\right)=$$

$$\displaystyle\left(\begin{array}[]{c c c}{\mbox{\boldmath$c$}}^{\,1}&\ldots&{{\mbox{\boldmath$c$}}^{\,N}}\end{array}\right)\left(\begin{array}[]{c c c }\delta_{1}{\mbox{\boldmath$\mu$}}&\ldots&\delta_{R+N_{c}}{\mbox{\boldmath$\mu$}}\end{array}\right).$$

$$\displaystyle\begin{cases}\frac{\partial{\color[rgb]{0,0,0}r}_{i}}{\partial x_{m}}\neq 0\ \ {\rm if}\ x_{m}\ {\rm influences\ reaction}\ i,\\ \frac{\partial{\color[rgb]{0,0,0}r}_{i}}{\partial x_{m}}=0\ \ {\rm otherwise.}\ \end{cases}{}$$

$$\displaystyle\nu=\left(\begin{array}[]{ccccc}1&-1&-1&1&0\\ 0&1&-1&1&-1\\ 0&0&1&-1&0\\ 0&0&1&-1&0\end{array}\right).$$

$$\displaystyle{\bf A}=\left(\begin{array}[]{cccc|cc}0&0&0&0&1&0\\ r_{2A}&0&0&0&1&0\\ r_{3A}&r_{3B}&0&0&0&1\\ 0&0&r_{4C}&r_{4D}&0&1\\ 0&r_{5B}&0&0&1&0\\ \hline\cr 0&0&1&-1&0&0\end{array}\right).$$

$${\delta }_{J^{\ast }}{x}_m$$

$${\delta }_{J^{\ast }}{r}_j$$

$$\lambda (\Gamma )\equiv -|\mathfrak{m}|+|\mathfrak{r}|-N(\mathfrak{r}).$$

$$\lambda (\Gamma )\equiv -|\mathfrak{m}|+|\mathfrak{r}|-N(\mathfrak{r})+N_c(\mathfrak{m}).$$

$$\{\mathrm{RasD},\mathrm{RasT},\mathrm{Raf},\mathrm{RafP},\mathrm{Mek},\mathrm{MekP},$$

$$\mathrm{MekPP},\mathrm{Erk},\mathrm{ErkP},\mathrm{ErkPP}\},$$

$$\displaystyle{\mbox{\boldmath$d$}}_{1}=(0,0,0,0,0,0,0,1,1,1)^{T}$$

$$\displaystyle{\mbox{\boldmath$d$}}_{2}=(0,0,0,0,1,1,1,0,0,0)^{T}$$

$$\displaystyle{\mbox{\boldmath$d$}}_{3}=(0,0,1,1,0,0,0,0,0,0)^{T}$$

$$\displaystyle{\mbox{\boldmath$d$}}_{4}=(1,1,0,0,0,0,0,0,0,0)^{T}$$

$$\displaystyle{\bf A}=\left(\begin{array}[]{cccccccccccc|cccc}r_{1,{\rm RasD}}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-1\\ 0&r_{2,{\rm RasT}}&0&0&0&0&0&0&0&0&0&0&0&0&0&-1\\ 0&r_{3,{\rm RasT}}&\pagecolor{LLLG}r_{3,{\rm Raf}}&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}-1&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}r_{4,{\rm RafP}}&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}-1&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}r_{5,{\rm RafP}}&\pagecolor{LLG}r_{5,{\rm Mek}}&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}-1&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}r_{6,{\rm MekP}}&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}-1&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}r_{7,{\rm RafP}}&\pagecolor{LLG}0&\pagecolor{LLG}r_{7,{\rm MekP}}&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}-1&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}r_{8,{\rm MekPP}}&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}-1&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}r_{9,{\rm MekPP}}&\pagecolor{LG}r_{9,{\rm Erk}}&\pagecolor{LG}0&\pagecolor{LG}0&\pagecolor{LG}0&\pagecolor{LG}-1&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LG}0&\pagecolor{LG}r_{10,{\rm ErkP}}&\pagecolor{LG}0&\pagecolor{LG}0&\pagecolor{LG}-1&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}r_{11,{\rm MekPP}}&\pagecolor{LG}0&\pagecolor{LG}r_{11,{\rm ErkP}}&\pagecolor{LG}0&\pagecolor{LG}-1&\pagecolor{LG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LG}0&\pagecolor{LG}0&\pagecolor{LG}r_{12,{\rm ErkPP}}&\pagecolor{LG}-1&\pagecolor{LG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ \hline\cr 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LG}-1&\pagecolor{LG}-1&\pagecolor{LG}-1&\pagecolor{LG}0&\pagecolor{LG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLG}-1&\pagecolor{LLG}-1&\pagecolor{LLG}-1&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0\pagecolor{LLG}&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLG}0&\pagecolor{LLLG}0&0\\ 0&0&\pagecolor{LLLG}-1&\pagecolor{LLLG}-1&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&\pagecolor{LLLG}0&0\\ -1&-1&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ \end{array}\right).$$

$$\displaystyle\begin{array}[]{c|cccccccccccc|cccc}&{k}_{1}&{k}_{2}&{k}_{3}&{k}_{4}&{k}_{5}&{k}_{6}&{k}_{7}&{k}_{8}&{k}_{9}&{k}_{10}&{k}_{11}&{k}_{12}&{l}_{1}&{l}_{2}&{l}_{3}&l_{4}\\ \hline\cr{\rm RasD}&-&+&0&0&0&0&0&0&0&0&0&0&0&0&0&+\\ {\rm RasT}&+&-&0&0&0&0&0&0&0&0&0&0&0&0&0&+\\ {\rm Raf}&-&+&-&+&0&0&0&0&0&0&0&0&0&0&+&-\\ {\rm RafP}&+&-&+&-&0&0&0&0&0&0&0&0&0&0&+&+\\ {\rm Mek}&-&+&-&+&-&+&-&+&0&0&0&0&0&+&-&-\\ {\rm MekP}&\pm&\pm&\pm&\pm&+&-&-&+&0&0&0&0&0&+&\pm&\pm\\ {\rm MekPP}&+&-&+&-&+&-&+&-&0&0&0&0&0&+&+&+\\ {\rm Erk}&-&+&-&+&-&+&-&+&-&+&-&+&+&-&-&-\\ {\rm ErkP}&\pm&\pm&\pm&\pm&\pm&\pm&\pm&\pm&+&-&-&+&+&\pm&\pm&\pm\\ {\rm ErkPP}&+&-&+&-&+&-&+&-&+&-&+&-&+&+&+&+\\ {}\hfil&&&&&&&&&&&&&&&&\\ {r}_{1}&+&+&0&0&0&0&0&0&0&0&0&0&0&0&0&+\\ {r}_{2}&+&+&0&0&0&0&0&0&0&0&0&0&0&0&0&+\\ {r}_{3}&+&-&+&+&0&0&0&0&0&0&0&0&0&0&+&+\\ {r}_{4}&+&-&+&+&0&0&0&0&0&0&0&0&0&0&+&+\\ {r}_{5}&\pm&\pm&\pm&\pm&+&+&-&+&0&0&0&0&0&+&\pm&\pm\\ {r}_{6}&\pm&\pm&\pm&\pm&+&+&-&+&0&0&0&0&0&+&\pm&\pm\\ {r}_{7}&+&-&+&-&+&-&+&+&0&0&0&0&0&+&+&+\\ {r}_{8}&+&-&+&-&+&-&+&+&0&0&0&0&0&+&+&+\\ {r}_{9}&\pm&\pm&\pm&\pm&\pm&\pm&\pm&\pm&+&+&-&+&+&\pm&\pm&\pm\\ {r}_{10}&\pm&\pm&\pm&\pm&\pm&\pm&\pm&\pm&+&+&-&+&+&\pm&\pm&\pm\\ {r}_{11}&+&-&+&-&+&-&+&-&+&-&+&+&+&+&+&+\\ {r}_{12}&+&-&+&-&+&-&+&-&+&-&+&+&+&+&+&+\\ \end{array},$$

$${\Gamma }_1=(\{\mathrm{Erk},\mathrm{ErkP},\mathrm{ErkPP}\},\{9,10,11,12\}).$$

$${\Gamma }_2$$

$${\Gamma }_3$$

Full text

PACS numbers

82.20.-w, 87.10.Ed, 87.18.-h

I Introduction

In living cells, there are many chemical reactions, and they forms complex networks such as metabolic networks and signal transduction networks. The dynamics of concentration of chemicals resulting from these networks is considered to be the origin of physiological functions. Today, a huge information about reaction networks is available in databases such as KEGG KEGG , Reactome Reactome and BioCyc BioCyc .

Dynamics and steady state of reaction networks are determined by various factors in systems; reaction-specific enzymes, input signals from outside of systems, and initial conditions in some cases. One of the standard approaches to elucidate the dynamics is a sensitivity experiment; for example, in the research of metabolism, changes in concentrations of metabolites induced by perturbation in amounts/activities of enzymes are examined Ishii (see Fig 1).

In previous theoretical studies Mochizuki ; monomolecular , under an implicit assumption that the system has no conserved concentrations, one of the authors and his collaborator showed that the sensitivity of steady state to reaction rate parameters, which correspond to enzyme activities/amounts, is determined only from the structure of the chemical reaction network (structural sensitivity analysis). Such a structural approach has a great advantage because it is almost impossible to measure the precise kinetics and parameters of chemical reactions in living cells. We then proved a novel theorem, called the law of localization OM . The law of localization characterizes subnetworks by non-positive indices (see $\lambda(\Gamma)$ in (60)) and identifies those which confine the effect of perturbations inside them. We call such subnetworks buffering structures.

One of the main purposes of this paper is to extend our previous method to reaction systems where some of concentrations are conserved during the dynamics. For example, in the MAPK pathway, while the ratio between activated and inactivated kinase concentrations changes after stimuli, the total number of them is conserved in time. In the presence of such conserved quantities, steady state concentrations and fluxes are influenced not only by reaction rate parameters but also by the initial values of conserved quantities. We take into account constraints coming from conserved quantities in the sensitivity analysis.

Another purpose of this paper is to explore biological meanings of buffering structures. Buffering structures provide the system with robustness to enzymatic fluctuations since they confine the effects inside them. Thus, we expect that possessing buffering structures are advantageous for living systems and networks with more buffering structures are selected during evolutionary process. In order to refine these expectations, we quantify robustness of reaction networks and compare robustness of Escherichia coli metabolism with artificial random networks.

The rest of the paper is organized as follows: In Sec. II, we generalize the method of structural sensitivity analysis. In Sec. III, we illustrate the method in a simple reaction network. In Sec. IV, we extend the law of localization to any reaction networks. In Sec. V, as applications, we study two signal transduction networks. The first network has conserved quantities and also include regulations from non-substrate chemicals. The second network also has conserved quantities. In Sec. VI, we propose network characteristics for robustness, and compare the E. coli network and random networks. The detailed explanation of the formulation and the proof of the law of localization are written in Appendix.

II Structural Sensitivity Analysis

The structural sensitivity analysis is a systematic method of determining sensitivity of steady states to rate parameter perturbations from reaction network information alone Mochizuki ; monomolecular . We generalize the structural sensitivity analysis to networks with conserved concentrations. In the generalized formulation, not only the sensitivity to rate parameter perturbations but also that to initial conditions on conserved concentrations are determined from network information.

We label chemical species by $m\,(m=1,\ldots,M)$ and reactions by $j\,(j=1,\ldots,R)$. In general, a macroscopic state of a spatially homogeneous chemical reaction system is specified by the concentrations $x_{m}(t)$ and obeys the following differential equations Mochizuki ; monomolecular ; F1 ; F2 ; OM

[TABLE]

Here, the $M\times R$ matrix $\nu$ is called the stoichiometric matrix: If the stoichiometry of the reaction $j$ among molecules $X_{m}$ is

[TABLE]

then the component $\nu_{mj}$ is defined as

[TABLE]

The reaction rate function $r_{j}$ is called a flux, which depends on the chemical concentrations $x$ and also on a reaction rate parameter $k_{j}$. We do not assume specific forms for the flux functions except that each flux is an increasing function of its substrate concentration;

[TABLE]

Below, we abbreviate $\frac{\partial r_{j}}{\partial x_{m}}$ evaluated at steady state as $r_{jm}$. Usual kinetics, such as the mass-action and the Michaelis-Menten kinetics, satisfies this condition.

In the context of metabolic reaction systems, $x_{m}$ is the concentration of the $m$-th metabolite, and the $j$-th reaction rate parameter $k_{j}$ corresponds the enzyme activity/amount catalyzing the $j$-th reaction.

We introduce notations about the kernel (right-null) and cokernel (left-null) spaces of $\nu$. We choose bases of the kernel and cokernel spaces, ${\rm ker}\ \nu$ and ${\rm coker}\ \nu$, as \{{\mbox{\boldmathc}}^{\alpha}\}_{(\alpha=1,\ldots,N)}and \{{\mbox{\boldmathd}}_{a}\}_{(a=1,\ldots,N_{c})}, where $N$ and $N_{c}$ denote their dimensions. The kernel space of $\nu$ corresponds to steady state fluxes Palsson ; FBA1 ; FBA2 ; FBA3 . Namely, steady state fluxes satisfy

[TABLE]

where $c^{\alpha}_{j}$ is the $j$-th component of the vector {\mbox{\boldmathc}}^{\alpha}, and $\mu^{\alpha}$ are coefficients. On the other hand, the cokernel space is related with conserved quantities. For every $a\in\{1,2,\ldots,N_{c}\}$,

[TABLE]

is constant in time, where $x$ is the vector $(x_{1},\ldots,x_{M})$. $l_{a}$ is a constant associated with {\mbox{\boldmathd}}_{a}.

In the previous papers Mochizuki ; monomolecular ; OM , the condition $N_{c}=0$ is assumed, and so steady state concentrations and fluxes are functions of rate parameters $k_{i}$. However, when $N_{c}>0$, steady state also depends on initial conditions on conserved concentrations, i.e. $\{l_{a}\}$. Therefore, in this case, there are two types of perturbations; the perturbation of the rate parameter, $k_{j^{*}}\rightarrow k_{j^{*}}+\delta k_{j^{*}}$ and that of the ${a^{*}}$-th conserved quantity, $l_{a^{*}}\rightarrow l_{a^{*}}+\delta l_{a^{*}}$, where $j={j^{*}}$ and $a={a^{*}}$ refer to the perturbed rate parameter and conserved quantity respectively.

To treat two types of perturbations in a unified way, we introduce generalized parameters $K_{J}\,(J=1,\ldots,R+N_{c})$ as

[TABLE]

We determine the concentration changes $\delta_{J^{*}}x_{m}$ and the flux changes $\delta_{J^{*}}r_{j}$ at steady state under the ${J^{*}}$-th perturbation, $K_{J^{*}}\rightarrow K_{J^{*}}+\delta K_{J^{*}}$ $({J^{*}}=1,\ldots,R+N_{c})$.

As shown in Appendix, responses to all perturbations are obtained simultaneously from the following matrix equation;

[TABLE]

where the vertical and horizontal lines indicate the structure of block matrices. The quantities above are defined as follows: First, ${\bf E}_{R}$, ${\bf E^{\prime}}_{N_{c}}$ are $R\times R$ and $N_{c}\times N_{c}$ diagonal matrices respectively, where the ${j^{*}}$-th component of ${\bf E}_{R}$ is given by $\frac{\partial r_{j^{*}}}{\partial k_{j^{*}}}\delta k_{j^{*}}$, and the ${a^{*}}$-th component of ${\bf E}_{N_{c}}$ is given by $\delta l_{a^{*}}$. For a fixed ${J^{*}}$, \delta_{J^{*}}\mbox{\boldmathx} and \delta_{J^{*}}\mbox{\boldmath\mu} in (12) are column vectors,

[TABLE]

where $\delta_{J^{*}}\mu^{\alpha}$ is the change of the coefficients in Eq. (5) under perturbation $K_{J^{*}}\rightarrow K_{J^{*}}+\delta K_{J^{*}}$. Finally, the matrix ${\bf A}$ is defined as

[TABLE]

Note that the matrix $\bf A$ is proved to be square.

By assuming that the matrix $\bf A$ is regular, Eq. (12) uniquely determines the sensitivity of chemicals, $\delta_{J^{*}}x_{m}$, and of the flux coefficients, $\delta_{J^{*}}\mu^{\alpha}$, as

[TABLE]

By using Eq. (5) and noting \mbox{\boldmathc}^{\alpha} is constant, $\delta_{J^{*}}\mu^{\alpha}$ determines the flux responses as

[TABLE]

or, in matrix notation,

[TABLE]

Note that \delta_{J^{*}}{\mbox{\boldmathr}} and {\mbox{\boldmathc}}^{\alpha} are $R$-dimensional column vectors, and \delta_{J^{*}}{\mbox{\boldmath\mu}} are $N$-dimensional column vectors.

Comments are in order. Firstly, practically, we are often interested in qualitative responses, (increased/decreased/invariant). For such discussions, assuming “overexpressions” $\delta K_{J}>0$, we can replace ${\bf E}_{R}$ and ${\bf E}^{\prime}_{N_{c}}$ by identity matrices. Therefore, we call ${\bf S}\equiv-{\bf A}^{-1}$ sensitivity matrix. Secondly, as a slight generalization, we can include nontrivial regulations such as allosteric effects by relaxing Eq. (4) as

[TABLE]

Such regulations add additional nonzero $r_{im}$ in the $\bf A$-matrix, but the response is still determined through Eq. (30). Finally, although we are using the terminology “perturbation”, the change of parameters $\delta K_{J^{*}}$ is not necessarily small for qualitative discussion of the sensitivity as long as the steady state persists for finite perturbations.

III Example network

We illustrate our method in a simple example shown in FIG. 2, which has $N_{c}=1$. See Appendix and OM for examples with $N_{c}=0$.

The stoichiometric matrix is given by

[TABLE]

$\nu$ has a cokernel vector {\mbox{\boldmathd}}_{1}=(0,0,1,-1) because the difference l_{1}\equiv{\mbox{\boldmathd}}_{1}\cdot{\mbox{\boldmathx}}=x_{C}-x_{D} is conserved in any reactions. $\bf A$ is given by

[TABLE]

From (12), the sensitivity is then calculated as

[TABLE]

where $R_{1}\equiv r_{4C}+r_{4D}$ and $R_{2}\equiv r_{2A}r_{3B}+r_{3A}r_{5B}$. The $J^{*}$-th column is the response to perturbations of $J^{*}$-th parameter and the vertical lines separate between perturbations of the rate parameters and the conserved quantity. For example, from the fifth column, we can see that the increase of the initial value of $l_{1}$ changes only $C$ and $D$, but does not affect either $A$, $B$ or all fluxes.

IV Law of localization

As proved in OM , for networks without conserved quantities, i.e. $N_{c}=0$, patterns of nonzero responses can be understood from network topology by using a theorem called the law of localization. The theorem is also useful for elucidating relevant pathways by combining with experimental measurements. In this section, we generalize the theorem into networks with conserved quantities, i.e. $N_{c}>0$.

First, we review the theorem for networks with $N_{c}=0$. For a given network, we consider a pair $\Gamma=({\mathfrak{m}},{\mathfrak{r}})$ of a chemical subset ${\mathfrak{m}}$ and a reaction subset ${\mathfrak{r}}$ satisfying the condition that ${\mathfrak{r}}$ includes all reactions influenced by metabolites in ${\mathfrak{m}}$ [see the condition Eq. (2’)]. We call $\Gamma=({\mathfrak{m}},{\mathfrak{r}})$ satisfying this condition output-complete. For a subnetwork $\Gamma=({\mathfrak{m}},{\mathfrak{r}})$, we define an index,

[TABLE]

Here, $|{\mathfrak{m}}|$ is the number of elements in ${\mathfrak{m}}$, $|{\mathfrak{r}}|$ the number of elements in ${\mathfrak{r}}$, and $N({\mathfrak{r}})$ the number of independent stoichiometric cycles in ${\mathfrak{r}}$. By a stoichiometric cycle in ${\mathfrak{r}}$ we mean any flux vector $c$ which satisfies the flux balance \nu\,{\mbox{\boldmathc}}=0 and nonzero components only within $\mathfrak{r}$. We call an output-complete subnetwork $\Gamma$ with $\lambda(\Gamma)=0$ as a buffering structure.

The law of localization states that the chemical concentrations and reaction fluxes outside of a buffering structure $\Gamma$ does not change under any rate parameter perturbations in $\mathfrak{r}$. In other words, all effects of perturbations of $k_{j^{*}}$ in $\mathfrak{r}$ are indeed localized within $\Gamma$. See Appendix for an illustration of the theorem for an example network with $N_{c}=0$.

Now, we state the theorem in the case of $N_{c}\geq 0$ (see Appendix for the proof). We replace the definition (60) of the index $\lambda$ by

[TABLE]

The additional contribution, $N_{c}({\mathfrak{m}})$, is the number of independent conserved quantities including at least one element in ${\mathfrak{m}}$. Note that the independencies of cycles and conserved quantities are defined in the vector spaces associated with $\mathfrak{r}$ and ${\mathfrak{m}}$ respectively (see Appendix B for a precise formal definition).

The generalized law of localization then states that the chemical concentrations and reaction fluxes outside of a buffering structure $\Gamma$ do not change under either perturbations of rate parameters or of conserved quantities in $\Gamma$.

We illustrate the generalized law of localization in the previous example network, shown in FIG. 2, in which there is a conserved quantity $l_{1}\equiv x_{C}-x_{D}$. $N_{c}({\mathfrak{m}})=1$ if we choose ${\mathfrak{m}}$ including $C$ or $D$, such as ${\mathfrak{m}}=\{C\},\{D\}$, $\{C,D\}$, or $\{A,C,D\}$. For the output-complete subgraph $\Gamma_{1}=(\{C,D\},\{4\})$, the index is $\lambda(\Gamma_{1})=-2+1-0+1=0$. For $\Gamma_{2}=(\{C,D\},\{3,4\})$, which has one cycle consisting of reactions $3,\ 4$, the index again vanishes; $\lambda(\Gamma_{2})=-2+2-1+1=0$. Therefore, these subnetworks are buffering structures. This is consistent for the result that the perturbation of the initial value of $x_{C}-x_{D}$ neither influences the concentrations of $A$ nor $B$. $\Gamma_{3}=(\{A,C,D\},\{2,3,4\})$ is another buffering structure with $\lambda(\Gamma_{3})=-3+3-1+1$, which explains why $x_{B}$ is insensitive to the perturbations of $k_{2},\ k_{3},\ k_{4}$, and $l_{1}$.

V Applications to biological networks

We apply the structural sensitivity analysis and the law of localization to two signal transduction pathways which have $N_{c}>0$.

V.1 Signal transduction 1: MAPK

We first consider the signal transduction network shown in FIG. 3.

The stoichiometric matrix $\nu$ of this system has a four-dimensional cokernel space corresponding to the total amounts of Ras, Raf, Mek, and Erk. Phosphorylated chemicals in the upper layer positively regulate phosphorylations in the lower one. Mathematically, this means that, for example, the arguments of the flux function $r_{7}$ is not the form $r_{7}(k_{7},{\rm MekP})$, but $r_{7}(k_{7},{\rm MekP,RafP})$, which additionally adds the component $r_{7,\rm RafP}$ in the $\bf A$ matrix (see Eq. (2’)).

To construct the matrix $\bf A$, we order the chemicals as

[TABLE]

where RasD/RasT denotes the bound state of Ras and ADP/ATP, and choose the following basis for cokernel vectors;

[TABLE]

The matrix $\bf A$ is given by

[TABLE]

Here, the gradations of color in ${\bf A}$ show block matrices corresponding to four buffering structures (see (106) and (107) below). The signs of the components of the sensitivity matrix ${\bf S}=-{\bf A}^{-1}$ are determined as

[TABLE]

where $+$, $-$ represent qualitative responses under perturbations associated with column indices, $k_{1},\ldots,l_{4}$. The symbol $\pm$ means that the sign depends on quantitative values of $r_{im}$.

The zero entries in (105) can be easily understood from the law of localization. The four square boxes in FIG. 3 indicate four buffering structures, forming a nested structure. The smallest buffering structure is

[TABLE]

This subnetwork has one conserved quantity, $l_{4}=x_{\rm Erk}+x_{\rm ErkP}+x_{\rm ErkPP}$, and so $N_{c}(\mathfrak{m})=1$. Therefore $\lambda(\Gamma_{1})=-3+4-2+1=0$. The law of localization then states that the perturbations of $k_{8},k_{9},k_{10},k_{11}$ and $l_{1}$ does not change the other part of the system, which explains the zeros appearing in the columns associated with $k_{9,10,11,12}$ and $l_{1}$ in Eq. (105).

The remaining buffering structures are

[TABLE]

These buffering structures explain the zero entries, and in particular, the nest of them explains the stair-like nonzero pattern in Eq. (105).

The nest of buffering structures implies that perturbations to an upper layer influence the lower layers of the signal transduction pathway.

V.2 Signal transduction 2: MAPK

We next study the signal transduction network shown in FIG. 4. This network was studied in sontag , where the regulation between the bottom two layers in FIG. 3 was modeled in detail as bound state formation in FIG. 4. They studied the sensitivity by using a clever manipulation of equations under the assumption of the mass-action kinetics. Here, by using structural sensitivity analysis, we derive the same result without assuming specific kinetics, which illustrates generality and usefulness of our method. We remark that, although this example is superficially similar to the one in FIG. 3, surprisingly, the result turns out to be more complex than the previous example.

This system has three conserved quantities associated with the total amounts of Erk, PTP, and Mek. We order the chemicals as

[TABLE]

and choose the basis of the cokernel space as

[TABLE]

The $\bf A$ matrix is given by

[TABLE]

The signs of responses are

[TABLE]

Note that this network has only one trivial buffering structure, i.e. the whole network. Accordingly, there are no vanishing entries in (139).

Following sontag , we focuse on total active Mek and Erk concentrations,

[TABLE]

and examine their responses to the perturbations of $l_{1},l_{2}$, which correspond to the total amounts of Erk and PTP with any form. In our method, these responses can be obtained by summing the associated rows in Eq. (139), which leads to

[TABLE]

Here we omit the common positive proportional constant, $({\rm Det}\,{\bf A})^{-1}>0$. We thus obtain the following qualitative responses,

[TABLE]

which agrees with the result obtained in sontag . The authors of sontag called the result (142) “paradoxical results”: While Eq. (142a) suggests that $x$ activates $y$, Eq. (142b) suggests that $x$ inhibits $y$.

We emphasize that while the argument of sontag is based on the mass-action type kinetics, we obtained the same conclusion valid for general (monotonically increasing) flux functions. Also, our systematic approach determines all responses simultaneously. We summarize the responses of $x,\,y$ for all perturbations;

[TABLE]

VI Characteristics for robustness based on network structures

In our previous study, we examined the carbon metabolism pathway of E. coli Mochizuki ; OM shown in FIG. 5, which is a major part of energy acquisition process. We found that the network has 17 buffering structures shown in FIG. 6 (see also Appendix for the list), and, in particular, that some of buffering structures coincide with subnetworks associated with biological functions such as the TCA cycle and the pentose phosphate pathway. These observations suggest that biological networks are selected in evolution and include many buffering structures, which provide robustness to enzymatic perturbations. In this section, in order to support this expectation, we compare robustness property between the E. coli network and artificial random networks.

Firstly, we introduce network characteristics that quantify robustness for reaction systems. One natural definition is the number of buffering structures, $N_{BS}$, which is more precisely defined as the number of buffering structures consisting of different sets of chemicals; we distinguish two buffering structures if they have at least one different metabolite. Note that we identify two buffering structures that have different sets of reactions even if they have the same metabolite set. Another one is the fraction of metabolites that exhibit zero responses under a randomly chosen enzymatic perturbation ($k_{j^{*}}\rightarrow k_{j^{*}}+\delta k_{j^{*}}$);

[TABLE]

By definition $0\leq{\mathscr{R}}\leq 1$. Networks with larger ${\mathscr{R}}$ and larger $N_{BS}$ are more robust. We emphasize that ${\mathscr{R}}$ and $N_{BS}$ are purely structural characteristics determined from network topology alone.

In order to generate random networks, we randomly choose $p\times R=48p$ reactions of the E. coli network and reconnect them to randomly chosen nodes of metabolites (see FIG. 7) alon . Here, $p$ ($0<p\leq 1$) is a fraction of rewired reactions. This is implemented by randomly choosing a column of the $M\times R$ stoichiometry matrix $\nu$, randomly reordering the $M$ components in the column, and repeat this procedure $48p$ times. Rewired networks are sometimes singular, i.e. ${\rm det}\,{\bf A}=0$, which cannot admit steady state, and sometimes unconnected, which is not suitable to compare with the E. coli network. We discard such networks. In this way, we constructed an ensemble of 3000 regular and connected networks for each value of $p$.

For these ensembles, we examined the robustness, ${\mathscr{R}}$ and $N_{BS}$, and the mean distance $D$, which is one of the most widely used network characteristics. Here, for a given network, the distance between two metabolites is defined as the smallest number of reactions that connects the two, and $D$ is the average value of the distances. More rigorously, $D$ here is defined for undirected networks where nodes represent metabolites and edges between two metabolites are drawn if the two metabolites are involved with the same reaction; for example, the distances between substrates and products of a reaction are defined as one.

The results are as follows. Firstly, FIG. 8 (a) is the distributions of ${\mathscr{R}}$ and $N_{BS}$ for $p=2/48,\,9/48,\,1$. The filled red circle represents the E. coli network, ${\mathscr{R}}_{E.\,coli}\simeq 0.57$ and N_{BS}\bigr{|}_{E.\,coli}\simeq 17. We can see a strong correlation between ${\mathscr{R}}$ and $N_{BS}$. Thus, the robustness based on ${\mathscr{R}}$ and $N_{BS}$ agree with each other statistically.

Secondly, FIG. 8 (b) shows the distribution of ${\mathscr{R}}$ and $D$ of a network. The blue circle represents the E. coli network, ${\mathscr{R}}_{E.\,coli}\simeq 0.57$ and $D_{E.\,coli}\simeq 2.45$. We can see that the distribution of ${\mathscr{R}}$ with fixed $D$ is widely spread for any $p$. This means that there are no significant correlations between ${\mathscr{R}}$ and $D$.

Finally, FIG. 8 (c) is the ensemble average of ${\mathscr{R}}$ and $D$. As we randomize the E.coli network by increasing $p$ from [math] to $1$, the robustness $\langle{\mathscr{R}}\rangle$ and the mean distance $\langle D\rangle$ tend to decrease monotonically. Here, $\langle\cdot\rangle$ denotes the average over the network ensemble with fixed $p$. We note that this tendency of $\langle{\mathscr{R}}\rangle$ and $\langle D\rangle$ does not necessarily mean that there is a correlation between ${\mathscr{R}}$ and $D$. In fact, we did not find such a correlation in FIG. 8 (b). We also confirmed these behaviors in a model like the Watts-Strogatz model WS (see Appendix).

One of the most remarkable results in FIG. 8 (a) is the peak of ${\mathscr{R}}$ at $p=0$, corresponding to the unrewired E. coli network. This peak means that even a small fraction of rewiring (for example $p=2/48$) lowers the robustness ${\mathscr{R}}$ drastically (from ${\mathscr{R}}_{E.\,coli}\simeq 0.57$ to $\langle{\mathscr{R}}\rangle\simeq 0.37$). The extraordinary robustness of the E. coli network suggests that the special topology, which is characterized by buffering structures, might be formed and selected under evolutionary pressures on the robustness.

While the robustness $\langle{\mathscr{R}}\rangle$ has a steep peak at $p=0$, $\langle D\rangle$ changes smoothly around $p=0$, as expected. In fact, other network quantities such as centrality, and degree distributions also change smoothly around $p=0$. The uncorrelated distributions shown in FIG. 8 (b) and the sharp peak of ${\mathscr{R}}$ at $p=0$ imply that ${\mathscr{R}}$ and $N_{BS}$ are completely novel characteristics for robustness, which cannot be captured by any conventional graph theoretical quantities.

VII Conclusions

In this paper, we generalized our previous formalism of structural sensitivity analysis and the law of localization into reaction systems with conserved quantities. Our generalized method can be applied into any biochemical systems, including signal transduction networks, metabolic systems, and protein synthesis, if the systems admit steady states.

We applied our method into two signal transduction networks with conserved quantities. While the authors in sontag studied the second network by assuming mass-action kinetics, we obtained the same conclusion as theirs without assuming specific kinetics such as mass-action types and the Michaelis-Menten kinetics. This illustrates how powerful and general our method is.

Our structural approach is also practically useful in experimental biology. In spite of the progress in biosciences, it is difficult to experimentally determine kinetics of biochemical reactions in living cells. Our method overcomes the difficulty because we determine qualitative sensitivity (increased/decreased/invariant) to perturbations only from network structures. By making use of this advantage, we can testify database information on networks systematically (see OM for more detailed discussions).

Finally, we investigated biological meanings of buffering structures by comparing E. coli network with random networks. We introduced two network characteristics measuring robustness of networks; the number of buffering structures $N_{BS}$ and the fraction of zero responses to perturbations $\mathscr{R}$. Based on them, we measured robustness of reactions systems. We found that even a partial rewiring deteriorates the robustness of the E. coli network drastically. E. coli network has special features that cannot be captured by other indexes than ours and realize extraordinarily robust compared with random networks. Our result suggests that the topology of the E. coli network might be selected under evolutionary pressures on robustness.

The proposed quantities for robustness, $N_{BS}$ and $\mathscr{R}$, are completely novel network characteristics because they are not correlated with conventional network characteristics, such as mean distance or degree distributions (see also Appendix). We will study more analytical aspects about the relation between robustness ($N_{BS}$ and $\mathscr{R}$) and network topology and hope to report on them in the near future.

This work was supported partly by the CREST, Japan Science and Technology Agency, and by iTHES Project/iTHEMS Program RIKEN, by Grant-in-Aid for Scientific Research on Innovative Area, “Logics of Plant Development,” Grant No. 25113001. We express our sincere thanks to Testuo Hatsuda, Michio Hiroshima, Yoh Iwasa, Masaki Matsumoto, Keiichi Nakayama, and Yasushi Sako for their helpful discussions and comments.

Appendix A Derivation of Equation (30)

As in many mathematical studies in metabolism like flux balance analysis FBA1 ; FBA2 ; FBA3 , we are focusing on steady state, which is characterized by

[TABLE]

Here, $\bar{x}_{m}$ represents the concentration of metabolite $m$ at steady state, which generally depends on (a subset of) the parameters $\{k_{j}\}$ and initial condition specified by $\{l_{a}\}$. We determine the sensitivity of each concentration and flux to the perturbations $k_{j^{*}}\rightarrow k_{j^{*}}+\delta k_{j^{*}}$ and $l_{a}\rightarrow l_{a}+\delta l_{a}$.

First, we determine the responses to perturbations of rate parameters. We choose one reaction $j=j^{*}$ in the system and perturb the parameter as $k_{j^{*}}\rightarrow k_{j^{*}}+\delta k_{j^{*}}$. The system goes to a new steady state, characterized by

[TABLE]

Here, $\delta_{jj^{*}}$ is the Kronecker delta symbol. Below we determine the concentration sensitivity $\delta_{j^{*}}x_{m}$ and the flux sensitivity $\delta_{j^{*}}r_{j}\equiv r_{j}(k_{j}+\delta_{jj^{*}}\delta k_{j^{*}},\bar{x}+\delta_{j^{*}}{x})-r_{j}(k_{j},\bar{x})$.

The above two equations imply that the flux change also satisfies $\sum_{j}\nu_{mj}\delta_{j^{*}}r_{j}=0$. Therefore, the vector $\delta_{j^{*}}{r}\in{\mathbb{R}}^{R}$ can be expanded in terms of a basis \{{\mbox{\boldmathc}}^{\alpha}\}_{\alpha=1,\ldots,N} of the kernel of $\nu$, where {{\mbox{\boldmathc}}^{\alpha}}\in{\mathbb{R}}^{R} and $N\equiv{\rm dim\,ker}\nu$;

[TABLE]

Here, ${c}^{\alpha}_{j}$ is the $j$-th component of the kernel vector {{\mbox{\boldmathc}}^{\alpha}}. Thus, the problem of determining the fluxes is equivalent to that of determining the coefficients $\delta_{j^{*}}\mu_{\alpha}$ of the kernel vectors.

As we commented in the main text, for the purpose of determining qualitative responses, we can assume that the perturbations are small. In the limit of infinitesimal perturbations $\delta k_{j^{*}}$, Taylor expansion of $\delta_{j^{*}}r_{j}$ yields

[TABLE]

where we abbreviate r_{jm}\equiv\frac{\partial r_{j}}{\partial x_{m}}\bigr{|}_{x={\bar{x}}}. Comparing Eqs. (150) and (151), we obtain

[TABLE]

When the cokernel vectors of $\nu$ exist, i.e. ${\rm dim\ coker}\,\nu\equiv N_{c}>0$ , Eq. (152) is not enough to determine the sensitivity, and we need additional constraints on the concentration changes. Let the constant vectors \{{\mbox{\boldmathd}}^{\,a}\}_{a=1,\ldots,N_{c}} be a basis of the cokernel space, where {\mbox{\boldmathd}}^{a}\in{\mathbb{R}}^{R}, and $N_{c}$ is the dimension of the cokernel space. Then the linear combinations

[TABLE]

where ({\mbox{\boldmathd}}_{a})_{m} denotes the $m$-th component of {\mbox{\boldmathd}}_{a}, are conserved in the dynamics of Eq. (1). This implies that the perturbed steady state depends on $\{l_{a}\}$. In order to make the problem of the sensitivity well-defined, we assume the perturbed system starts with the same initial condition as the unperturbed system. Then the concentration changes need to satisfy

[TABLE]

for all $a=1,\ldots,N_{c}$.

Eqs. (152) and (154) determine the response to the perturbation ${j^{*}}$. In matrix notation, these can be written as

[TABLE]

where the matrix $\bf A$ is defined in Eq. (25), and {\mbox{\boldmathe}}_{j^{*}}\equiv(0,\ldots,\underset{{j^{*}}{\rm{\mathchar 45\relax}th}}{\underbrace{\frac{\partial r_{j^{*}}}{\partial k_{j^{*}}}\biggr{|}_{x=\bar{x}}\delta k_{j^{*}}}},\ldots,0)^{T}\in{\mathbb{R}}^{R}.

We note that the matrix $\bf A$ is square, namely the identity $M+N=R+N_{c}$ holds. This follows from the well-known identity for the Fredholm index, ${\rm dim\,ker}\,\nu-{\rm dim\,coker}\,\nu=R-M$ for any $M\times R$ matrix $\nu:{\mathbb{R}}^{R}\rightarrow{\mathbb{R}}^{M}$.

Next, we discuss responses to perturbation of conserved quantities, or initial concentrations. We choose a particular conserved quantity, $a=a^{*}$, and consider the perturbation $l_{a^{*}}\rightarrow l_{a^{*}}+\delta l_{a^{*}}$. Eqs. (150), (151), (152), and (154) are replaced by

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

From Eq. (166) and (167), we obtain the matrix equation,

[TABLE]

where $\bf A$ is the same as that in Eq. (25), and the column vector {\mbox{\boldmathe}}_{a^{*}} is defined as {\mbox{\boldmathe}}_{a^{*}}\equiv(0,\ldots,\underset{{a^{*}}{\rm{\mathchar 45\relax}th}}{\underbrace{\delta l_{a^{*}}}},\ldots,0)^{T}\in{\mathbb{R}}^{N_{c}}.

If we write the results of the perturbations for $j^{*}=1,\ldots,R$ and ${a^{*}}=1,\ldots,N_{c}$, given by Eq. (163) and (176), we obtain Eq. (30) in the main text.

Appendix B Proof of the law of localization

Firstly, we write the precise definitions of $N({\mathfrak{r}})$ and $N_{c}({\mathfrak{m}})$ appearing in Eq. (61). For a chemical subset ${\mathfrak{m}}$ and a reaction subset ${\mathfrak{r}}$, we can respectively associate the following vector spaces $V({\mathfrak{r}})$ and $V_{c}({\mathfrak{m}})$,

[TABLE]

Here, $P^{\mathfrak{r}}$ is an $R\times R$ projection matrix onto the space associated with ${\mathfrak{r}}$ defined as

[TABLE]

In other words, $V({\mathfrak{r}})$ are vectors {\mbox{\boldmathv}}\in{\mathbb{R}}^{R} with component support in $\mathfrak{r}$. Similarly, $P^{\mathfrak{m}}$ is an $M\times M$ projection matrix on the space associated with ${\mathfrak{m}}$. Then, we define $N({\mathfrak{r}})$ and $N_{c}({\mathfrak{m}})$ as the dimensions of these vector spaces;

[TABLE]

The intuitive meaning for this definition is explained in the main text. Note that $N({\mathfrak{r}})={\rm dim\,ker}(\nu P^{\mathfrak{r}})-R+|\mathfrak{r}|$ and $N_{c}({\mathfrak{m}})={\rm dim\,coker}\,\nu-{\rm dim\,coker}(P^{\bar{\mathfrak{m}}}\nu)+|\mathfrak{m}|$, where $P^{\bar{\mathfrak{m}}}\equiv\hat{\bf 1}_{M\times M}-P^{\mathfrak{m}}$ is a projection matrix on the space associated with the complementary chemicals $\bar{\mathfrak{m}}$ of ${\mathfrak{m}}$.

Now we prove the theorem. Suppose that $\Gamma=({\mathfrak{m}},{\mathfrak{r}})$ is a buffering structure, namely an output-complete subnetwork satisfying $\lambda(\Gamma)=0$. As discussed below, by choosing appropriate bases of the kernel and the cokernel of $\nu$ and arranging the orders of the column and row indices of the matrix ${\bf A}$, we can always rewrite $\bf A$ into the form,

[TABLE]

Since we are assuming ${\rm Det}\,{\bf A}\neq 0$, the upper-left block in Eq. (187) is generally vertically long or square; i.e. $\lambda(\Gamma)\geq 0$. The condition $\lambda(\Gamma)=0$ means that it is square.

The structure of block matrices in Eq. (187) can be obtained by collecting the indices associated with $\Gamma=({\mathfrak{m}},{\mathfrak{r}})$ into the upper-left corner: The column indices at the upper left block consist of the chemicals in ${\mathfrak{m}}$ followed by the basis vectors of the kernel space $V({\mathfrak{r}})$ associated with $\mathfrak{r}$. The row indices consist of the reactions in $\mathfrak{r}$ followed by the basis vectors of the cokernel space $V_{c}({\mathfrak{m}})$ associated with ${\mathfrak{m}}$. Thus, from this construction, the upper-left block has the size $(|{\mathfrak{r}}|+N_{c}({\mathfrak{m}}))\times(|{\mathfrak{m}}|+N({\mathfrak{r}}))$. We need to prove that the lower-left block vanishes completely. First, all $r_{jm}$ with $j\notin{{\mathfrak{r}}}$ and $m\in{\mathfrak{m}}$, which would appear in the lower-left block, vanish by the assumption that the subnetwork $\Gamma$ is output-complete; reaction rates $r_{j}(k_{j},x)$ with $j\notin{\mathfrak{r}}$ do not depend on the concentration of chemicals in $\mathfrak{m}$. It remains to show that {\mbox{\boldmathc}}^{\alpha}\in V({\mathfrak{r}}) and {\mbox{\boldmathd}}_{a}\in V_{c}(\mathfrak{m}) do not have nonzero entries in the lower-left block. But this directly follows from their definitions, Eq. (177). This completes the proof of the structure given in Eq. (187).

After arranging the matrix $\bf A$ as in Eq. (187), it is easy to prove the law of localization. We first prove the theorem for the reaction rate perturbations, i.e. $\delta_{j^{*}}{x}_{m}=0$ and $\delta_{j^{*}}r_{j^{\prime}}=0$ for ${j^{*}}\in{{\mathfrak{r}}}\not\ni j^{\prime}$, and $m\notin{{\mathfrak{m}}}$. As explained in the body of the paper, the concentration change $\delta_{j^{*}}{x}_{m}$ is proportional to $({\bf A}^{-1})_{m{j^{*}}}\propto\ {\rm Det}\,{\hat{\bf A}}^{{j^{*}},m}$, where ${\hat{\bf A}}^{{j^{*}},m}$ is the minor matrix obtained by removing reaction row ${j^{*}}$ and chemical column $m$ from the matrix $\bf A$. Noting that ${j^{*}}\in{{\mathfrak{r}}}$ belongs to the raw indices of the upper part of $\bf A$, and $m\notin{{\mathfrak{m}}}$ belongs to the column indices of of the right part of $\bf A$, ${\rm Det}\,{\hat{\bf A}}^{{j^{*}},m}=0$ holds for ${j^{*}}\in{{\mathfrak{r}}}$ and $m\notin{{\mathfrak{m}}}$ because the nonzero block at the upper left of the minor ${\hat{\bf A}}^{{j^{*}},m}$, which was originally square in Eq. (187), now becomes horizontally long. This proves $\delta_{j^{*}}{x}_{m}=0$ for ${j^{*}}\in{{\mathfrak{r}}}$ and $m\notin{{\mathfrak{m}}}$. It remains to show $\delta_{j^{*}}r_{j^{\prime}}=0$ for all ${j^{*}}\in{{\mathfrak{r}}}\not\ni j^{\prime}$. Noting $r_{j^{\prime}}$ depends on the outside chemicals $x_{m}$ with $m\notin\mathfrak{m}$ because of the output-completeness of $\Gamma=({\mathfrak{m}},{\mathfrak{r}})$, Eq. (151) becomes $\delta_{j^{*}}r_{j^{\prime}}=\sum_{m\notin{\mathfrak{m}}}\frac{\partial r_{j^{\prime}}}{\partial x_{m}}\delta_{j*}x_{m}$. Then, the chemical insensitivity $\delta_{j^{*}}x_{m}=0$ for $m\not\in\mathfrak{m}$ also means the flux insensitivity $\delta_{j^{*}}r_{j^{\prime}}=0$.

Similarly, for the perturbations of conserved quantities, we can prove the chemical insensitivity, $\delta_{a^{*}}x_{m}=0$ for all $m\notin\Gamma$ and the flux insensitivity $\delta_{a^{*}}r_{j^{\prime}}=0$ for all $j^{\prime}\notin{{\mathfrak{r}}}$, under any perturbation of conserved quantities in $\Gamma$, that is, any perturbation of l_{a^{*}}={\mbox{\boldmathd}}_{a^{*}}\cdot{\mbox{\boldmathx}} for {\mbox{\boldmathd}}_{a^{*}}\in V_{c}({\mathfrak{m}}). From (176), $\delta_{a^{*}}x_{m}$ is proportional to $({\bf A}^{-1})_{ma^{*}}\propto\ {\rm Det}\,{\hat{\bf A}}^{a^{*},m}$, the determinant of the minor matrix obtained by removing the row of the $a^{*}$-th conserved quantity. Noting that $m\notin\Gamma$ belongs to the column indices of the right part of $\bf A$ and the conserved quantity $l_{a^{*}}$ associated with {\mbox{\boldmathd}}_{a^{*}}\in V_{c}({\mathfrak{m}}) belongs to the row indices of the upper part of $\bf A$, we can show ${\rm Det}\,{\hat{\bf A}}^{a^{*},m}=0$ and $\delta_{a^{*}}x_{m}=0$, which leads to the flux insensitivity $\delta_{a^{*}}r_{j^{\prime}}=0$, as in the above argument for the reaction rate perturbations.

This proves the law of localization. $\square$

Appendix C Example network

We illustrate the computation of the sensitivity analysis for the network consisting of $R=6$ reactions and $M=4$ chemicals, shown in FIG. 9.

The stoichiometric matrix $\nu$ is

[TABLE]

Noting that $N_{c}=0$, the matrices A and $\bf S$ are

[TABLE]

[TABLE]

Then, from Eqs. (30) and (36), the responses of chemical concentrations and fluxes are

[TABLE]

and

[TABLE]

We can see that only the perturbation to the input rate, corresponding to the 1st column in Eq. (206), affects all chemicals and fluxes. The perturbations to reactions $2,3,5$ only decrease the concentrations of the substrates $A,B,D$ respectively. The perturbation of reaction $4$ decreases the concentrations $D,A,B$ along the cycle downward of the perturbation (see FIG. 9, and the 4th column of $\bf S$). The perturbation of reaction 6 does not change the further downstream but change $A,B,C,D$ in the cycle.

The law of localization can be applied as follows. Some of the buffering structures are shown in FIG. 10. The smallest buffering structures are $\Gamma_{1}=(\{A\},\{2\}),\ \Gamma_{2}=(\{B\},\{3\}),\ \Gamma_{3}=(\{D\},\{5\})$, which all sassily $\lambda(\Gamma_{i})=-1+1-0=0$. In addition, the network has two larger ones, $\Gamma_{4}=(\{A,B,D\},\{2,3,4,5\})$ (with $\lambda(\Gamma_{4})=-3+4-1=0$), $\Gamma_{5}=(\{A,B,C,D\},\{2,3,4,5,6\})$ (with $\lambda(\Gamma_{5})=-4+5-1=0$). $\Gamma_{4}$ is the minimum buffering structure including reaction 4. Then, the law of localization predicts that the nonzero response to perturbation of reaction 4 should be limited within $\Gamma_{4}$, which is observed in the 4th column in Eq. (206). Similarly, the response to perturbation of reaction 6 is explained by $\Gamma_{5}$.

Appendix D E. coli central metabolism

D.0.1 List of reactions

1: Glucose + PEP $\rightarrow$ G6P + PYR.

2: G6P $\leftarrow$ F6P.

3: F6P $\rightarrow$ G6P.

4: F6P $\rightarrow$ F1,6P.

5: F1,6P $\rightarrow$ G3P + DHAP.

6: DHAP $\rightarrow$ G3P.

7: G3P $\rightarrow$ 3PG.

8: 3PG $\rightarrow$ PEP.

9: PEP $\rightarrow$ 3PG.

10: PEP $\rightarrow$ PYR.

11: PYR $\rightarrow$ PEP.

12: PYR $\rightarrow$ AcCoA + CO2.

13: G6P $\rightarrow$ 6PG.

14: 6PG $\rightarrow$ Ru5P + CO2.

15: Ru5P $\rightarrow$ X5P.

16: Ru5P $\rightarrow$ R5P.

17: X5P + R5P $\rightarrow$ G3P + S7P.

18: G3P + S7P $\rightarrow$ X5P + R5P.

19: G3P + S7P $\rightarrow$ F6P + E4P.

20: F6P + E4P $\rightarrow$ G3P + S7P.

21: X5P + E4P $\rightarrow$ F6P + G3P.

22: F6P + G3P $\rightarrow$ X5P + E4P.

23: AcCoA + $\rightarrow$ CIT.

24: CIT $\rightarrow$ ICT.

25: ICT $\rightarrow$ 2${\rm\mathchar 45\relax}$KG + CO2.

26: 2-KG $\rightarrow$ SUC + CO2.

27: SUC $\rightarrow$ FUM.

28: FUM $\rightarrow$ MAL.

29: MAL $\rightarrow$ OAA.

30: OAA $\rightarrow$ MAL.

31: PEP + CO2 $\rightarrow$ OAA.

32: OAA $\rightarrow$ PEP + CO2.

33: MAL $\rightarrow$ PYR + CO2.

34: ICT $\rightarrow$ SUC + Glyoxylate.

35: Glyoxylate + AcCoA $\rightarrow$ MAL.

36: 6PG $\rightarrow$ G3P + PYR.

37: AcCoA $\rightarrow$ Acetate.

38: PYR $\rightarrow$ Lactate.

39: AcCoA $\rightarrow$ Ethanol.

40: R5P $\rightarrow$ (output).

41: OAA $\rightarrow$ (output).

42: CO2 $\rightarrow$ (output).

43: (input) $\rightarrow$ Glucose.

44: Acetate $\rightarrow$ (output).

45: Lactate $\rightarrow$ (output).

46: Ethanol $\rightarrow$ (output).

D.0.2 List of buffering structures

The E. coli network exhibits the following 17 different buffering structures $\Gamma_{i}=({\mathfrak{m}}_{i},{\mathfrak{r}}_{i})$ ($i=1,\ldots,17$).

$\Gamma_{1}=(\{\rm Glucose\},\{1\})$,

$\Gamma_{2}=(\{\rm Glucose,PEP,G6P,F6P,F1,6P,DHAP,G3P,3PG,\\ PYR,6PG,Ru5P,X5P,R5P,S7P,E4P,AcCoA,OAA,CIT,\\ ICT,2{\rm\mathchar 45\relax KG},SUC,FUM,MAL,CO2,Glyoxylate,Acetate,\\ Lactate,Ethanol\},\{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,\\ 16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,\\ 34,35,36,37,38,39,40,41,42,44,45,46\})$,

$\Gamma_{3}=(\{\rm F1,6P\},\{5\})$,

$\Gamma_{4}=(\{\rm DHAP\},\{6\})$,

$\Gamma_{5}=(\{\rm G3P,X5P,S7P,E4P\},\{7,17,18,19,20,21,22\})$,

$\Gamma_{6}=(\{\rm 3PG\},\{8\})$,

$\Gamma_{7}=(\{\rm Glucose,PEP,3PG,PYR,AcCoA,OAA,CIT,ICT,\\ 2{\rm\mathchar 45\relax KG},SUC,FUM,MAL,CO2,Glyoxylate,Acetate,Lactate,\\ Ethanol\},\{1,8,9,10,11,12,23,24,25,26,27,28,29,30,31,\\ 32,33,34,35,37,38,39,41,42,44,45,46\})$,

$\Gamma_{8}=(\{\rm X5P,S7P,E4P\},\{17,18,19,20,21\})$ ,

$\Gamma_{9}=(\{\rm CIT\},\{24\})$,

$\Gamma_{10}=(\{\rm 2{\rm\mathchar 45\relax KG}\},\{26\})$,

$\Gamma_{11}=(\{\rm SUC\},\{27\})$ ,

$\Gamma_{12}=(\{\rm FUM\},\{28\})$

$\Gamma_{13}=(\{\rm Glyoxylate\},\{35\})$,

$\Gamma_{14}=(\{\rm X5P,R5P,S7P,E4P\},\{17,18,19,20,21,40\})$ ,

$\Gamma_{15}=(\{\rm Acetate\},\{44\})$,

$\Gamma_{16}=(\{\rm Lactate\},\{45\})$,

$\Gamma_{17}=(\{\rm Ethanol\},\{46\})$.

Appendix E The degrees and the values of robustness in rewired E. coli networks

We define the degree of the $m$-th metabolite as the number of reactions with which the $m$-th metabolite participates is involved as a substrate or product. Here, we computed the variance of degrees and $\mathscr{R}$ for each rewired network, and investigated whether there exist any correlation between these two quantities. We note that the average of degrees in a network are the same for all rewired networks because the rewiring procedure preserves the total number of reactions. FIG. 11 shows the distributions of the variance of degree and $\mathscr{R}$ for 556 rewired networks when $p=4$. We did not observe any strong correlation between them, as we mentioned in the main text.

Appendix F The Watts-Strogatz-like model

Here, we consider a model of random reaction networks similar to the Watts-Strogatz model WS .

The model consists of $M$ chemicals and reactions along a circular clockwise pathway and randomly directed reactions inside the circle. The end points of internal reactions are randomly rewired from the next neighbors to other chemicals with probability $p$. There is an inflow and an outflow on the circle. In fact, ${\mathscr{R}}$ and $N_{BS}$ are independent of the positions of the inflow and the outflow. We note that, unlike the original Watts-Strogatz model, the reactions along the circle are not rewired, which guarantees the existence of the inverse of $\bf A$ in Eq. (11).

We generated 1000 random networks for each ensemble labeled by $p$. Two figures in FIG. 13 show the distributions of ${\mathscr{R}},N_{BS}$ and of ${\mathscr{R}},D$. We can see that while $N_{BS}$ and ${\mathscr{R}}$ are all positively correlated for any value of $p$, ${\mathscr{R}}$ and $D$ are not correlated significantly.

FIG. 14 shows the result of $\langle{\mathscr{R}}\rangle$ and $\langle D\rangle$ for various $p$. As we randomize the network by increasing $p$, the mean distance $\langle D\rangle$ and $\langle{\mathscr{R}}\rangle$ decrease monotonically. In contrast to the E. coli network, there is no strong peak of the robustness ${\mathscr{R}}$ around $p=0$ in the Watts-Strogatz-like model.