A structural approach to understanding enzymatic regulation of chemical reaction networks

Mochizuki, Atsushi

doi:10.1042/BCJ20210545

In living cells, chemical reactions are connected by sharing their products and substrates, and form complex systems, i.e. chemical reaction network. One of the largest missions in modern biology is to understand behaviors of such systems logically based on information of network structures. However, there are series of obstacles to study dynamical behaviors of complex network systems in biology. For example, network structure does not provide sufficient information to determine details of the dynamical behaviors. In this review, I will introduce a novel mathematical theory, structural sensitivity analysis, by which the responses of reaction systems upon the changes in enzyme activities/amounts are determined from network structure alone. The patterns of responses exhibit characteristic features, localization and hierarchy, depending on the topology of the network. The theory also shows that ranges of enzymatic regulations are governed by a mathematical law characterized by local topology of substructures. These findings imply that the network topology is one of the origins of biological robustness.

Introduction

The biological functions of a cell arise from a large set of chemical reactions. In the living cell the reactions define an inter-connected large system, where products of one reaction act as substrates of other reactions. The relationship between chemicals and reactions are represented by chemical reaction networks, where a node represents a species of chemicals and a reaction arrow represents a state-transition of chemicals.

It is widely believed that the dynamics of chemicals based on such complex networks are the origin of biological functions. One of the largest missions in modern biology is to obtain logical understanding of the behaviors of systems based on experimentally identified networks. On the other hand, in contrast with the rapid increase in information of network structures, we have a limited understanding of dynamical behaviors of network systems. There are a variety of obstacles impeding attempts to study the behaviors of biological systems based on the reaction networks.

One of the difficulties is the observation of dynamic processes. It is still difficult to observe the dynamics of the chemicals with sufficient time resolution. Most of the data obtained by present experimental methods are snapshots of molecular concentrations rather than time tracks. The second problem is the reliability of the reaction network itself. At present the information of network structures are possibly incomplete in many studies of biological systems because of the complexity. The problem is fundamental because we can never exclude the possibility that unknown reactions may take an important role in the focal phenomena.

The third and largest problem is that the information on the network structure alone is not sufficient to determine the resulting dynamics. The reaction edges only provide qualitative information on state-transition between chemical species in the system. They lack essential quantitative details like the reaction functions, parameter values of reaction rates, and initial states. A typical approach to the difficulty is modeling: developing mathematical models for a given network under assumptions of functions, parameters and initial conditions, and trying to reproduce the observed behaviors. Such numerical simulations therefore rely on many unverified assumptions as to the reaction functions and their dozens or hundreds of unknown parameters. In general, numerical parameter identification does not seem to be a viable option given the lack of highly-sampled time series data collected in most biological experiments.

In this review, we introduce a new mathematical approach to determine the behavior of complex chemical reaction networks structurally. Although total and detailed dynamical behaviors cannot be determined solely from network structure, important aspects of dynamical behaviors can be determined from network structure alone. Our focal property here is the responses of the systems to the changes in enzyme activities in the system, i.e. how reaction systems are regulated by the modulation of enzymatic activities. By the method ‘structural sensitivity analysis', the qualitative change in the concentration of chemicals in the system is determined from information on the reaction network, only, without any assumptions on either the specific functions modeling the reactions, or their reaction parameters. We will show that (1) nonzero responses of the system are localized in finite regions in networks, and show characteristic patterns [1,2]. We also show that (2) a general law, which we call ‘law of localization’, directly governs the patterns of nonzero responses. Any substructure satisfying an equation: (#. chemical sp.) − (#. reactions) + (#. cycles) = 0, where the right-hand-side is an analog of Euler characteristic, is a buffering structure, i.e. any perturbation of reactions inside the substructure does not change concentrations and fluxes outside the structure [3,4].

Chemical reaction networks

The biological function of cells is fundamentally originated from the chemical reactions of biomolecules in a cell. Many types of chemical reactions work in the cell, and are connected to form a chained network. For example, Figure 1 shows the central metabolic pathway, which is a universally observed in various biological species from bacteria to humans. In this pathway, sugar (glucose) is decomposed into smaller organic substances, and chemical energy is extracted as ATP in the process. The rate of such biochemical reactions is controlled by reaction-specific organic catalysts known as enzymes. It is known that many chemical reactions working in a living cell are linked in a chain in a way that shares a substrate and a product to form a huge network. It is thought that the physiological functions of cells emerge from the dynamics of this entire system. In addition, it is thought that the physiological function is controlled via changes in the amount or activity of the enzymes.

Reaction network of the central metabolic pathway of bacteria (Escherichia coli).

Figure 1.

Includes 28 chemicals and 46 chemical reactions. Arrows indicate changes in substances due to chemical reactions. The numbers attached to the arrows indicate the chemical reaction, and the character strings at both ends of the arrow indicate the metabolites. The colors correspond to Figure 6. ([3], modified from [5]).

View large Download slide

Reaction network of the central metabolic pathway of bacteria (Escherichia coli).

Includes 28 chemicals and 46 chemical reactions. Arrows indicate changes in substances due to chemical reactions. The numbers attached to the arrows indicate the chemical reaction, and the character strings at both ends of the arrow indicate the metabolites. The colors correspond to Figure 6. ([3], modified from [5]).

Information on chemical reactions in living organisms is collected and organized in large-scale databases, and complex networks such as those shown in Figure 1 can be seen in various biological studies. These networks are naturally constructed in the process that the identified reactions are connected by sharing the substrate and the product. However, in contrast with the amount of information on reactions, we have surprisingly limited understanding about the dynamics of chemical reaction network. Also, our current collection of knowledge may not be incomplete.

One experimental approach to understand the behavior of reaction systems is sensitivity: perturbing the system experimentally and measuring the response. In particular, the enzymes catalyzing each reaction in the network are perturbed via knockdown or overexpression, and changes in the concentration of chemical substances are measured as the response. A method of measuring many metabolites in the system at once using mass spectrometry has been actively used since the 2000s, and is called metabolome analysis. These experiments are not just performed to measure responses to perturbations. They can also be understood as pseudo-control experiments in which changes in the activity of enzymes which control physiological functions in organisms are induced artificially.

Determining fluxes of reactions in a system is a little more technical compared with measuring the concentration of chemicals. Metabolic flux analysis is a method to estimate the fluxes of reactions from the changes of patterns of labeled chemicals [6]. One or more substrates labeled by radioisotope are introduced into cells, concentrations of labeled chemicals are measured after some time by mass spectrometry or nuclear magnetic resonance. From the distribution of labeled chemicals, fluxes of reactions in focal system are estimated.

On the other hand, interpretation of the results obtained from these perturbation experiments is extremely difficult as well. For example, there is a pioneering study in which various enzymes were knocked down to the central metabolic pathway of E. coli and metabolome analysis was performed for each [5]. In this study, in many knockdown experiments, only a few metabolites increased or decreased in concentration, and most metabolites showed no significant change. Even when changes are observed, the response patterns are controversial. Based on these results, Ishii et al. considered the possibility of the existence of unknown bypass reactions not described in the database. We need a criterion to discuss the results of perturbation experiments based on solid logic rather than intuition. In the following, we introduce a theoretical method to determine the response of a system to a change in the activity or quantity of an enzyme only from the network structure.

Dynamics of chemical reaction networks

The theory of the dynamics of chemical reaction systems has a long history, yet the fundamental idea is simple. Accordingly, the change in the concentration of each substance in the system over time is determined by the difference in inflow and outflow caused by reactions.

For example, in the reaction system shown in Figure 2a, the concentration of molecule A increases with the reaction rate of reaction 1 and decreases with the reaction rate of reaction 2. For each molecule, the change in concentration over time is related with the difference in the reaction rates by the differential equation system:

\frac{d u_{A}}{d t} = r_{1} - r_{2}

1a

\frac{d u_{B}}{d t} = r_{2} - r_{3}

1b

Here,

u_{m} \cdot (m = A, B)

is the concentration of the chemical in the system, and

r_{j} \cdot (j = 1, 2, 3)

is the reaction rate of each reaction. Eq. (1) can alternatively be expressed as follows:

\frac{d}{d t} (\begin{matrix} u_{A} \\ u_{B} \end{matrix}) = (\begin{matrix} 1 & - 1 & 0 \\ 0 & 1 & - 1 \end{matrix}) (\begin{matrix} r_{1} \\ r_{2} \\ r_{3} \end{matrix})

2

The dynamics of a chemical reaction system then can generally be described as follows:

\frac{d}{d t} (\begin{matrix} u_{1} \\ ⋮ \\ u_{M} \end{matrix}) = (\begin{matrix} ν_{11} & \dots & ν_{1 J} \\ ⋮ & ⋮ \\ ν_{M 1} & \dots & ν_{M J} \end{matrix}) (\begin{matrix} r_{1} \\ ⋮ \\ r_{J} \end{matrix})

3

Here,

(u_{1}, \dots, u_{M})^{T}

are the concentration vectors of chemicals in the system, and

(r_{1}, \dots, r_{J})^{T}

are the reaction rate vectors.

ν = (ν_{m j})

is termed the stoichiometric matrix. Each row and column of the stoichiometric matrix corresponds to a substance and a reaction, respectively. Each element

ν_{m j}

of this matrix indicate increases or decreases of molecule m resulting from a reaction j occurring once.

ν

is determined from the structure of the reaction network. The reaction rate function

r_{j}

is commonly modeled using the mass reaction kinetics

r_{2} = k_{2} u_{A}

or the Hill function

r_{2} = k_{2} u_{A}^{n} / (K^{n} + u_{A}^{n})

⁠.

Sensitivity analysis of a small reaction network.

Figure 2.

(a) Example of a chemical reaction system. A and B represent the chemicals, whereas 1,2, and 3 represent the reactions. (b) Examples of perturbations and responses to chemical reaction systems. Red arrowheads indicate reactions with increasing reaction parameter. Red circles (dotted lines) indicate decrease in molecular concentrations. Black circles and black arrows indicate no change in the steady state in response to perturbation. (c) Dynamics of molecular concentrations and reaction rates. The red arrowhead indicates the time point at which the perturbation was applied.

View large Download slide

Sensitivity analysis of a small reaction network.

(a) Example of a chemical reaction system. A and B represent the chemicals, whereas 1,2, and 3 represent the reactions. (b) Examples of perturbations and responses to chemical reaction systems. Red arrowheads indicate reactions with increasing reaction parameter. Red circles (dotted lines) indicate decrease in molecular concentrations. Black circles and black arrows indicate no change in the steady state in response to perturbation. (c) Dynamics of molecular concentrations and reaction rates. The red arrowhead indicates the time point at which the perturbation was applied.

Although not discussed in detail here, there is an idea to facilitate mathematical analysis by assuming a form of mass action in the reaction rate function [7,8]. Some mathematical conclusions on dynamics based on models using mass action kinetics are summarized in [9]. For example, it is proved that the concentrations of chemicals remain nonnegative for nonnegative initial concentrations as long as the solution exists. Feinberg and his colleagues, proved that, in a strongly connected chemical reaction network, 0 value of an index called ‘Deficiency', which is defined solely from the structure of a network, is a sufficient condition for the existence of a unique steady state solution [10–13].

Structural sensitivity analysis

Let us focus on the steady state of the reaction system defined in Eq. (3). The assumption of a steady state is because in general, the rates of chemical reactions are much higher than those rates of expression of genes encoding enzymes. I introduce a method to determine the steady state response to the change in each enzyme from network structure [1–4]. Analyses to determine the response to a perturbation are called sensitivity analysis in general.

There are some studies to analyze sensitivity of chemical reaction networks. Kacser and Burns [14] and Heinrich and Rapoport [15] independently proposed a mathematical idea, which has been called ‘metabolic control analysis' later [16,17]. The analysis provides a mathematical framework to determine the sensitivity of a single pathway of chemical reactions and of some cases of branched system. The metabolic control analysis has been applied to some examples of small chemical reaction networks. I will introduce a different and simpler mathematical framework to study the sensitivity by a function-free approach, which enables to relate the structure of the network and the sensitivity responses of the system directly.

The method introduced here is called ‘structural sensitivity analysis' since it is based on network structure alone. Note that although the expression ‘perturbation of the enzyme' is used here, sensitivity to any constant that affects the reaction rate can be given by the same formula. We do not assume a specific functional form for the reaction rate function $r_{j}$ ⁠. We assume only two trivial facts: (1) Concentration dependence: the rate of each reaction depends on the concentration of the substrate or the regulating molecule, (2) Reaction parameters: each reaction rate has a constant that corresponds to the activity or amount of enzyme. For (1), the concentration of the substrate molecule is naturally contained in the arguments of the reaction rate function, but it is also known that molecules other than the substrate may regulate the enzyme activity, as it is called allosteric regulation. In the following, if allosteric regulation exists, it is explicitly annotated on the network structure. In both cases, the dependences are expressed by non-zero partial derivative $\partial r_{j} / \partial u_{m}$ of reaction rate $r_{j}$ with respect to the concentration of each molecule, $u_{m}$ ⁠.

Here, we present only the method and results of the structural sensitivity analysis (See Appendix A for the derivation of the method). The following two-step calculation allows us to determine any molecular concentration and reaction rate responses to the perturbation of any reaction parameters in the system at once. In the following, we consider the case where there is no conserved quantity in the system. This corresponds to the case where the stoichiometric matrix is full rank (⁠

rank (ν) = M

⁠) (Appendix A). First, a matrix is constructed from the chemical reaction network:

A := (\begin{matrix} \frac{\partial r_{1}}{\partial u_{1}} & \dots & \frac{\partial r_{1}}{\partial u_{M}} \\ ⋮ & ⋮ \\ \frac{\partial r_{J}}{\partial u_{1}} & \dots & \frac{\partial r_{J}}{\partial u_{M}} \end{matrix} | \begin{matrix} - c^{1} & \dots & - c^{N} \end{matrix}) .

4

The left part

\partial r_{j} / \partial u_{m}

of matrix

A

is the partial derivative of the reaction rate

r_{j}

with respect to the chemical concentrations

u_{m}

⁠, and the right part

c^{n}

is the basis vector of the null space of the stoichiometric matrix

ν

(⁠

ker (ν) = {c \in R^{J} | ν c = 0}

⁠). Each row of matrix

A

corresponds to each reaction

(j = 1, \dots, J)

and the first half of the column corresponds to each molecular species

(m = 1, \dots, M)

⁠. In the second half of the columns of matrix

A

⁠,

- c^{n}

is arranged in each column

(n = 1, \dots, N)

⁠. Note that N is the dimension of null space

v

⁠. (If there is no conserved quantity, then

J = M + N

⁠.) Each component

\partial r_{j} / \partial u_{m}

of the left part shows the dependence of each reaction on each molecule, and is zero where reaction j is independent of molecule m. If a reaction j depends on a molecule m (e.g. it is the substrate or the regulator of the reaction), it then obtains a non-zero value. Nonzero component of the left part is explicitly denoted below with

r_{j m} (= \partial r_{j} / \partial u_{m})

⁠. The right part

c^{n}

⁠, i.e. the basis vector in the null space of matrix

ν

⁠, corresponds to the flow path in the steady state determined from the network.

For example, for the reaction system shown in Figure 2a (Eq. (2)), stoichiometric matrix

ν

is full rank (⁠

rank (ν) = 2

⁠),

dim (ker (ν)) = 3 - rank (ν) = 1

⁠, and matrix

A

is given by:

A = (\begin{matrix} 0 & 0 & - 1 \\ r_{2 A} & 0 & - 1 \\ 0 & r_{3 B} & - 1 \end{matrix}) .

5

Next, the inverse of matrix

A

is calculated. Then, the change in the concentration of each molecule (or change in the rate of each reaction) when a perturbation is given to the parameter of each reaction is given as each component of the inverse matrix:

S \equiv (\begin{matrix} \begin{matrix} δ_{1} u_{1} & \dots & δ_{J} u_{1} \\ ⋮ & ⋮ \\ δ_{1} u_{M} & \dots & δ_{J} u_{M} \end{matrix} \\ \bar{\begin{matrix} δ_{1} μ^{1} & \dots & δ_{J} μ^{1} \\ ⋮ & ⋮ \\ δ_{1} μ^{N} & \dots & δ_{J} μ^{N} \end{matrix}} \end{matrix}) \propto - A^{- 1} .

6

Here

δ_{j} u_{m}

indicates the change in the concentration of molecule m upon perturbation of the parameter of reaction j.

δ_{j} μ^{n}

denotes the change in steady flow

c^{n}

upon perturbation of the parameter of reaction j. The change in the reaction rate vector

r

upon perturbation of the parameter of reaction j is also given as:

δ_{j} r = \sum_{n = 1}^{N} δ_{j} μ^{n} c^{n} .

7

Here, zero or nonzero of each component of sensitivity matrix

S

is determined from the distribution of zero (and nonzero) components of matrix

A

⁠. As stated above, the distribution of the zero components of matrix

A

is determined solely from the structure of the reaction network. Thus, the response of the system to perturbations of reaction constants, i.e. the presence or absence of changes in the concentration of each molecule or the rate of each reaction, can be determined from the network structure alone. Furthermore, by assuming that each reaction rate is a monotonically increasing or decreasing function of the concentration of each molecule, the sign of the change, i.e. increase or decrease in concentration or rate, can also be determined from the network structure alone in many cases.

For the dynamics defined in Eq. (2) (Matrix

A

in (5)), the sensitivity matrix

S

is given by:

S = (\begin{matrix} 1 / r_{2 A} & - 1 / r_{2 A} & 0 \\ 1 / r_{3 B} & 0 & - 1 / r_{3 B} \\ 1 & 0 & 0 \end{matrix}) .

8

Here, each column refers to each of perturbed reaction parameters (from 1 to 3), and row 1, 2, and 3 represent the response of A, B, and the steady flow, respectively.

Figure 2b illustrates the case when parameter $k_{2}$ of reaction 2 is increased. This perturbation corresponds to an experiment in which the activity or expression level of the enzyme catalyzing reaction 2 is increased. One might expect an increase in the concentration of the product molecule B and an increase in the rate of the perturbed reaction. However, the conclusions of the structural sensitivity analysis are different. As shown in column 2 of (8), the only response to the increase in the enzyme activity of reaction 2 is decrease in molecule A. Neither the concentration of B nor any of the three reaction rates at steady state change in response to the perturbation.

The theory does not depend on the shape of the reaction rate function, but when an increasing function of substrate is chosen for the reaction rates and the dynamics are calculated, the behavior shown in Figure 2b is observed. Immediately after increasing the enzyme activity, the reaction rate $r_{2}$ and the concentration $u_{B}$ of the downstream molecule increase, but this is only temporary. Reaction rate $r_{2}$ quickly decreases to its pre-perturbation level, whereas the concentration of substrate, $u_{A}$ decreases. At the steady state, only the concentration of substrate A is observed to change in response to perturbation. Qualitatively the same behavior is observed regardless of the form of reaction rate function $r_{j}$ (as long as it monotonically increases with respect to substrate concentration). We may interpret that the change in the reaction parameter is compensated by the change in the concentration of the substrate molecule, and the substrates and reaction rates downstream do not change at all.

Let us give another example: For the chemical reaction system shown in Figure 3, the dynamics of the concentrations of A, B, and C are given by the differential equation system:

\frac{d}{d t} (\begin{matrix} u_{A} \\ u_{B} \\ u_{C} \end{matrix}) = (\begin{matrix} 1 & - 1 & 1 & 0 & 0 \\ 0 & 1 & - 1 & - 1 & 0 \\ 0 & 0 & 0 & 1 & - 1 \end{matrix}) (\begin{matrix} r_{1} \\ r_{2} \\ r_{- 2} \\ r_{3} \\ r_{4} \end{matrix}) .

9

Note that stoichiometric matrix

ν

is full rank,

rank (ν) = 3

and

dim (ker (ν)) = 5 - rank (ν) = 2

so

\ker (ν)

has two basis vectors. Matrix

A

is given by:

A = (\begin{matrix} 0 & 0 & 0 & - 1 & 0 \\ r_{2 A} & 0 & 0 & - 1 & - 1 \\ 0 & r_{- 2 B} & 0 & 0 & - 1 \\ 0 & r_{3 B} & 0 & - 1 & 0 \\ 0 & 0 & r_{4 C} & - 1 & 0 \end{matrix}) .

10

Here, each row corresponds to each reaction (⁠

1, 2, - 2, 3, 4

⁠), and each column corresponds to chemical A, B, C or to basis vectors of

ker (ν)

⁠.

r_{j m}

’s are nonzero components of the partial derivative of the reaction rate

r_{j}

with respect to the chemical concentrations

u_{m}

⁠. For example,

r_{2 A}

in (2, 1) indicates that the rate of reaction 2 depends on the concentration of A. The two basis vectors correspond to the pathways that run through reactions

1, 2, 3, 4

from inflow to outflow (column 4) and to the pathways of the cycle in the system consisting of reactions 2 and

- 2

(column 5).

Example of chemical reaction system.

Figure 3.

(0) A–C denote three speciess of chemicals in the system, 1,2–2,3,4 denote reactions, and the arrows denote state changes due to these reactions. (1–4) Response patterns to the increase in each reaction parameter in the chemical reaction system shown in (0). The results of structural sensitivity analysis (11) are graphically represented. The red arrowhead shows a perturbed reaction. Red letters indicate changes in molecular concentration. Red upper and lower cases denote increase and decrease in concentration, respectively. Red arrows indicate changes in reaction rate. Red bold and dotted arrows denote increase and decrease in reaction rate, respectively. Black letters and black arrows indicate unchanged concentrations of molecules and reactions, respectively.

View large Download slide

Example of chemical reaction system.

(0) A–C denote three speciess of chemicals in the system, 1,2–2,3,4 denote reactions, and the arrows denote state changes due to these reactions. (1–4) Response patterns to the increase in each reaction parameter in the chemical reaction system shown in (0). The results of structural sensitivity analysis (11) are graphically represented. The red arrowhead shows a perturbed reaction. Red letters indicate changes in molecular concentration. Red upper and lower cases denote increase and decrease in concentration, respectively. Red arrows indicate changes in reaction rate. Red bold and dotted arrows denote increase and decrease in reaction rate, respectively. Black letters and black arrows indicate unchanged concentrations of molecules and reactions, respectively.

Sensitivity matrix

S

in this case is given by:

S = (\begin{matrix} \frac{r_{- 2 B} + r_{3 B}}{r_{2 A} r_{3 B}} & - \frac{1}{r_{2 A}} & \frac{1}{r_{2 A}} & - \frac{r_{- 2 B}}{r_{2 A} r_{3 B}} & 0 \\ \frac{1}{r_{3 B}} & 0 & 0 & - \frac{1}{r_{3 B}} & 0 \\ \frac{1}{r_{4 C}} & 0 & 0 & 0 & - \frac{1}{r_{4 C}} \\ 1 & 0 & 0 & 0 & 0 \\ \frac{r_{- 2 B}}{r_{3 B}} & 0 & 1 & - \frac{r_{- 2 B}}{r_{3 B}} & 0 \end{matrix}) .

11

The responses of the system to the increase in the parameter of each reaction are obtained using (11) as follows (Figure 3):

Reaction 1: The concentrations of A, B, and C increase. The fluxes of the pathways through the system and the cycle pathway consisting of reactions 2 and $- 2$ increase.

Reaction $2$ ⁠: The concentration of A decreases. The concentrations of the other molecules are unchanged. Fluxes of all reactions are unchanged.

Reaction $- 2$ ⁠: The concentration of A increases, whereas those of the other molecules remain unchanged. The fluxes of the cycle pathway consisting of reactions $2$ and $- 2$ increases.

Reaction 3: The concentrations of A and B decrease, whereas that of C does not change. The flux of the cycle pathway consisting of reactions $2$ and $- 2$ increases.

Reaction 4: The concentration of C decreases, whereas those of the other molecules remain unchanged. Fluxes of all reactions are unchanged.

Note that perturbation responses of reactions 2 and $- 2$ are asymmetrical, and that the effect of perturbation of reaction 3 is transmitted two steps upstream. Thus, characteristic response patterns are observed depending on the topology of the network and the perturbed reactions.

Application (1): plant metabolism

Below, we will introduce an example of applying structural sensitivity analysis to an actual biological system. Ferjani et al. found that their mutants had metabolic abnormalities [18]. In other words, the sucrose production is significantly reduced in the early development of the mutant plants. Sucrose is an important sugar molecule used for transport between tissues in plants. On the other hand, it was revealed that this mutant loses the function of H⁺-pyrophosphatase, an enzyme that degrades pyrophosphate (PPi). PPi is naturally formed by the hydrolysis of ATP into AMP in many biological processes in cells, and the accumulation of them to toxic levels disrupts several common biosynthetic pathways and is lethal. Figure 4a shows the pathway by which sucrose is synthesized from fructose in early plant development. Pyrophosphate is synthesized as a by-product of two reversible reactions on this pathway. Based on these findings, the following story was initially envisioned. In individuals in which the pyrophosphate degrading enzyme is disrupted, the concentration of pyrophosphate would increase, which would make slowing down of the forward direction of the two reversible reactions on the sucrose synthesis pathway. Therefore, the rate of sucrose synthesis downstream of the pathway may decrease. At first glance, this story sounds reasonable.

Structural sensitivity analysis of the plant metabolism system.

Figure 4.

(a) Metabolic pathway for sucrose synthesis from fructose in early plant development. Molecules that are supposed to be abundant in the cytoplasm, such as phosphate groups, are not described. The change in the concentration of each metabolite when the PPi-degrading enzyme was knocked down was determined by structural sensitivity analysis and shown by the color of each molecule [18]. Red: increase, blue: decrease, gray: no change. The dotted arrow indicates an example of modification of the newly added reaction. Adding this reaction does not change the qualitative response to PPi-degrading enzyme perturbations. (b) The effect on the sucrose response when various reactions are added to the plant metabolic pathway given in Figure 4a. Green: The addition of these reactions does not alter the sucrose concentration response to the disruption of PPi-degrading enzyme. An arrow with arrowheads at both ends means that the result does not depend on the reaction direction or whether it is reversible reaction or not. Blue: When any of these four reactions is added, the sucrose concentration decreases by the destruction of PPi-degrading enzyme. Red: When this reaction is added, the sucrose concentration rather increases with the destruction of PPi degrading enzyme.

View large Download slide

Structural sensitivity analysis of the plant metabolism system.

(a) Metabolic pathway for sucrose synthesis from fructose in early plant development. Molecules that are supposed to be abundant in the cytoplasm, such as phosphate groups, are not described. The change in the concentration of each metabolite when the PPi-degrading enzyme was knocked down was determined by structural sensitivity analysis and shown by the color of each molecule [18]. Red: increase, blue: decrease, gray: no change. The dotted arrow indicates an example of modification of the newly added reaction. Adding this reaction does not change the qualitative response to PPi-degrading enzyme perturbations. (b) The effect on the sucrose response when various reactions are added to the plant metabolic pathway given in Figure 4a. Green: The addition of these reactions does not alter the sucrose concentration response to the disruption of PPi-degrading enzyme. An arrow with arrowheads at both ends means that the result does not depend on the reaction direction or whether it is reversible reaction or not. Blue: When any of these four reactions is added, the sucrose concentration decreases by the destruction of PPi-degrading enzyme. Red: When this reaction is added, the sucrose concentration rather increases with the destruction of PPi degrading enzyme.

To verify this story, Ferjani et al. applied structural sensitivity analysis to the sucrose synthesis pathway in plants. Figure 4a shows the change in concentration (increase/decrease/no change) as a result of disruption of PPi-degrading enzymes (pyrophosphatase). Contrary to expectations, disruption of pyrophosphatase does not result in a decrease in sucrose. Certainly, the concentration of pyrophosphate increases by the destruction of pyrophosphatase. At the same time, the concentration of the four sugars upstream in the pathway increases, and the concentration of UDP-glucose, which is the product of the reaction releasing pyrophosphate, decreases. So far, the story is true. However, the concentration of sucrose-6P, product of complex formation of fructose-6P and UDP-glucose, does not change. The concentration of sucrose downstream of sucrose-6P is also unchanged. So far as using the reaction network shown in Figure 4a, this result is always true without depending on the choice of reaction rate function or reaction parameters. On the other hand, in actual plants, a decrease in sucrose concentration has certainly been observed in mutants of pyrophosphatase. This means that the network structure in Figure 4a must be reconsidered and revised.

The dotted arrow in Figure 4a is an example of a modification of the network. A network with the additional reaction was assumed, and the response to the destruction of pyrophosphatase was calculated using the network. In this case, the response of the concentration of each sugar is exactly the same as before the addition was applied. This modification cannot explain the observed phenomenon. There are many other possible modifications. Ferjani et al. added each of reactions to the network, and analyzed the response to pyrophosphatase disruption in these modified networks. These results are summarized in Figure 4b.

The majority of the added reactions do not alter the response of the system to pyrophosphatase disruption. Only the networks modified to include any of the four reactions shown in blue show the observed behavior, namely the decrease in sucrose to the destruction of pyrophosphatase. On the other hand, the addition of the reaction shown in red causes the opposite behavior to the observation, that is, sucrose concentration rather increases by the destruction of pyrophosphatase.

From these results, Ferjani et al. predict that one of the reactions shown in blue exists and the reaction shown in red does not exist in the actual plant metabolism system. This prediction was later confirmed experimentally. In other words, it was confirmed by experimental measurement that the concentration of the products of the reaction shown in blue changed. Thus, by combining the structural sensitivity analysis with the knowledge of the network and the observed behavior, it is possible to predict existence of an unknown reaction specifically on a network system.

Law of localization

In the previous sections, we showed that the qualitative response (increase, decrease, or no change) of the chemical reaction system to the enzyme perturbation is determined from the structure of the network alone. In each example, the effects of perturbations on the reaction parameters are confined in a finite range of the system in steady state, i.e. the influence of perturbations is restricted in some chemicals and some reactions on the networks. From these results, we can say that the common intuitive argument that ‘upstream perturbations affect downstream' is not correct.

This kind of behavior is indeed universally observed as previously reported in reaction networks [19] and regulatory networks [20]. That is, the effect of perturbations on the reaction parameter is often confined to a limited part of the reaction network. There is a clear boundary on the network between the range affected by perturbations and that unaffected. Okada and Mochizuki discovered the mathematical law behind this phenomenon that in fact governs the behavior of reaction systems [3,4].

Let us first select a substructure from the reaction network arbitrarily, i.e. we choose a subset of chemicals and a subset of reactions under the condition that (i) the reactions dependent on chemicals in the chemical subset must be included in the reaction subset. This condition is called ‘output complete'. After satisfying this condition, another reaction may be added to the subset. Note that an output complete substructure includes all of emanating arrows from the selected chemicals, but may not include the substrate (or regulator) chemicals of the selected reactions. The second condition is that (ii) the number of chemical species, the number of reaction types, and the number of cycles in the substructure satisfy the following equation: (number of chemical species) − (number of reactions) + (number of cycles) = 0. A substructure satisfying both (1) and (2) is then called ‘buffering structure'. The buffering structure presents the following properties:

Theorem (Law of Localization)

A perturbation given to a reaction parameter inside a buffering structure affect only the steady state of chemical concentrations and reaction rates inside the buffering structure. The steady state of chemical concentrations and reaction rates outside the buffering structure is completely unaffected by the perturbations inside the buffering structure.

The proof is obtained only from the distribution of the zero components of the matrix $A$ (See Appendix B for details).

The example in Figure 5a is the simplest buffering structure. The structure containing one molecule and one reaction emanating from the molecule satisfies the condition of the buffering structure (⁠ $1 - 1 + 0 = 0$ ⁠). In general, if there is only one reaction arrow emanating from a molecule, any perturbation applied to the reaction will affect only the substrate molecule. Figure 5b shows the case of the plant metabolism system. PPi is drawn in two places, but since they are the same molecule, the number of molecules inside the red frame is 6. Before the modification (without including the green reaction), the number of reactions is 11 as the reversible reaction counts one forward and one backward. There are five reversible reactions in the red frame conforming small cycles. The structure in the red frame satisfies the conditions of the buffering structure (⁠ $6 - 11 + 5 = 0$ ⁠). In other words, although we discussed the effect of perturbation to PPi-degrading enzymes in Section 5, not only that, but also the perturbation to any of the reaction parameters in the red frame in Figure 5b has no effect on the concentration of sucrose-6P or sucrose.

Examples of buffering structures.

Figure 5.

(a) Corresponds to Figure 2b and Figure 4a. The crosses indicate perturbations upon the parameters on the reactions. Note that in (b), the PPi drawn in two places are the same molecule. Without depending on the presence or absence of the reaction shown in green, the index value of the condition of buffering structure remains 0.

View large Download slide

Examples of buffering structures.

(a) Corresponds to Figure 2b and Figure 4a. The crosses indicate perturbations upon the parameters on the reactions. Note that in (b), the PPi drawn in two places are the same molecule. Without depending on the presence or absence of the reaction shown in green, the index value of the condition of buffering structure remains 0.

In Figure 4b, various network modifications were made, but most of them did not alter the sucrose response at all. This can be understood as follows. In Figure 5b, adding the reaction shown in green increases the number of reactions by one. However, at the same time, a new cycle structure is also created through the addition of the green reaction. Therefore, the value of the conditional expression remains 0 (⁠ $= 6 - 12 + 6$ ⁠), and the buffering structure does not change.

There were only five network modifications in which the perturbation of PPi-degrading enzyme affected the concentration of sucrose (Figure 4b). All of these were additions of reactions emanating outward from the molecules inside the buffering structure shown in Figure 5b. According to the condition of buffering structure, we see that only the addition of this kind of reaction is able to break the condition. In other words, the five reactions described in Figure 4b are all of possible modifications that alter the sucrose response, i.e. they are not only sufficient but also necessary conditions for altering sucrose response.

Note the fact that any substructure satisfies the conditions for a buffering structure alone creates a localized response to perturbations. Properties of no change to enzyme fluctuations, in other words robustness, may arise simply from the condition of numbers of objects in subnetworks without intentionally creating a system. Okada and Mochizuki created a network randomly to see how often buffering structures would be formed [4], and found that buffering structures appear frequently in random networks.

Application (2): central metabolic pathway

Okada and Mochizuki analyzed responses of central metabolic pathway of E. coli (Figure 1) to perturbations of reaction parameters systematically [3]. Nonzero responses of the concentrations of chemicals to the perturbations of each reaction parameter in the system are summarized in Figure 6. The figure shows the localized and hierarchical property in the response patterns. When, for example, any one of the six reactions in the yellow box is perturbed, only three molecules in the yellow box respond. Similarly, when reaction 7 (shown in red) is perturbed, four molecules in the red and yellow boxes respond; whereas when reaction 40 (shown in orange) is perturbed, four molecules in the orange and yellow boxes respond. When the reaction shown in blue is perturbed, there are 11 cases where only one molecule responds, and in other cases all blue molecules respond. When the reaction shown in green is perturbed, green, yellow, red, orange, and blue molecules respond. Thus, not only are the responses to perturbations of various reactions respectively localized, the response ranges show inclusion relations with each other. Hence, perturbation responses in a reaction network are localized and hierarchical.

The perturbation response network of the E. coli central metabolic pathway obtained from the analysis.

Figure 6.

The numbers at the bottom of each box indicate the reaction parameters that were perturbed. The chemical name in each box and in all boxes followed by an arrow extending downward indicates the chemicals whose concentration changes in response to a given perturbation. Colors correspond to those in Figure 1 [3].

View large Download slide

The perturbation response network of the E. coli central metabolic pathway obtained from the analysis.

The numbers at the bottom of each box indicate the reaction parameters that were perturbed. The chemical name in each box and in all boxes followed by an arrow extending downward indicates the chemicals whose concentration changes in response to a given perturbation. Colors correspond to those in Figure 1 [3].

The hierarchy of responses patterns to perturbations observed in the central metabolic pathway of E. coli can be understood as the nest of buffering structures. This network has 17 buffering structures nested between one another. In Figure 6, a buffering structure corresponds to chemicals and reactions shown in each box and in the boxes followed by the downward arrows.

The nest of buffering structures may have functional and evolutionary significances. The function of the buffering structure is analogically understood with a fire wall. It is a mechanism that keeps the influence of fluctuations occurring inside from being transmitted outside. In other words, the central metabolic pathway of E. coli has 17 nested firewall-like structures that absorb given fluctuations inside and prevent transmission of them to the outside. The characteristic property of the system creates strong homeostasis or robustness outside against changes in enzyme amount and activity inside.

The following hypotheses may be possible to understand the biological significance of the mechanism to create external homeostasis in the central metabolic pathway in E. coli. The central metabolic pathway is a circuit in which chemical energy is extracted from nutrients (sugars) taken from the outside of a cell, and the activities or amounts of enzymes in the circuit have to be regulated always according to the nutritional environment. On the other hand, there are multiple circuits that are indispensable for maintaining cell functions, such as amino acid synthesis pathways and lipid synthesis pathways, connecting to the central metabolic pathway. Physiological problems will arise, if the regulation of enzyme activity in the central metabolic pathway influences the pathways for maintaining cell functions. For this reason, buffering structures may exist in the central metabolic pathway. The whole chemical reaction network contains multiple pathways with different physiological functions, which are inter-connected via reactions. Many other buffering structures may thus be used in various parts of chemical reaction networks to achieve functional independence (modularity) between connected pathways.

Discussion

In this review, I introduced a new mathematical analysis for a chemical reaction network, structural sensitivity analysis. By the analysis, responses of chemical concentrations and reaction rates in a reaction system to a perturbation of reaction parameter is determined from the network structure alone. I also introduced the principle of law of localization and the concept of buffering structure that has special significance for the range of influence of parameter change on the network.

We expect that the structural sensitivity analysis will be used as a tool to promote the researches of exchanging theoretical predictions and experimental verifications in life sciences. For example, we may calculate the sensitivity of a reaction system by the structural sensitivity analysis based on a network information on a database. At the same time, we may measure the response of the system by giving perturbations upon enzymes in the system experimentally. If there exists any disagreement between theoretical predictions and experimental measurements, it directly indicates the existence of inconsistencies between the database and the real network. As the theory is model-free, we do not have to consider the possibility of inappropriate choices of reaction functions or parameter values.

After finding inconsistencies between the database and the real network, the theory is also used to find candidates of unknown reactions or regulations in the networks. For example, we may calculate the sensitivity for hypothetical modifications which introduce, or omit, some reactions or regulations, and may then compare with experimental results. If such modifications considerably enhance agreement with experimental results, then we may expect the modified network to actually represent the real systems. The theory may be useful also to guide experimental design for small molecule treatment. For example, we may select potential drug targets in the metabolic system by the theory. The candidate targets may be further experimentally validated and lead to potential treatments.

A buffering structure is a substructure that absorbs effects of inside fluctuations and prevent these influences from being transmitted to the outside under the steady state condition. The law of localization may denote a fundamental principle for creating homeostasis and adaptation in living systems. One of the characteristics behaviors of life is called adaptation: a system may temporarily respond to changes in the environment such as light and nutrition, but return to the original state eventually. Most of earlier studies on adaptation have assumed that negative feedback acts from downstream to upstream in the signal pathway to understood the transient response. The law of localization implies that even without this kind of active regulation, adaptation can still occur if some part of the network structure satisfies certain topological conditions.

The law of localization also implies that the extent of enzymatic control is determined solely by the local structure of the network. This may then be used as a guideline to identify huge systems in step-by-step manner from a small substructure to the whole networks. We search for a candidate buffering structure based on the information in the database. Next, we measure the response of the system by giving a perturbation to this candidate structure experimentally. If, contrary to the expectation, the perturbation influences the outside of the candidate buffering structure, it indicates that some inconsistency should exist inside (or immediately near) the candidate structure, because the condition of buffering structure depends only on the local structure. Then we search for unknown reactions or regulations inside this structure by some experimental methods, and update the network structure using the newly-obtained information. By repeating this process of prediction and verification multiple times, the correct network structure can be determined in a stepwise manner from a small structure to a large structure.

Similarly, depending on the aim of study, it may also be possible to use the network structure as a guide to determine the scope of systems of interest in advance. For example, if the control of biological functions by enzymes is being studied, the focus should be on the smallest buffering structure that contains the enzymatic reaction of interest, since if the network structure is correct, it is mathematically guaranteed that the effects of fluctuations will not transmitted outside the identified buffering structure.

The recent progresses in the structural theory show that the buffering structures govern not only the sensitivity but also bifurcation property of the reaction systems [21,22]. Discontinuous transition of behavior of a system induced by a continuous change in parameters are called bifurcation, and have been studied well in the field of dynamical systems theory. Analytical study of bifurcation behaviors of a large system is usually very difficult, and numerical analysis needs a large computational load. Okada et al. showed that bifurcation behaviors in a complex network can be studied by decomposing it into smaller subnetworks based on a topological condition on networks [21,22]. For each subnetwork, the condition for bifurcation occurrence and chemicals exhibiting bifurcation behaviors are determined on the network. In other words, a steady state bifurcation may take place modularly in a reaction system, and the extent of the bifurcation behaviors are determined by the buffering structure.

Since the methodology presented here is based on only network information without any assumptions on the model, it may have a strong application potential to assist experimental biology studies. A large effort is paid on research to understand dynamics of network systems in life sciences. Using structure sensitivity analysis, we will be able to understand the properties of living systems by linking them to the network structure.

Appendix A: Mathematics of structural sensitivity analysis

In a chemical reaction system, the time variation of the concentration of a substance is generally expressed in form of a differential equations system as follows:

\frac{d u_{m}}{d t} = \sum_{j = 1}^{J} ν_{m j} r_{j} (k_{j}; u) . (A .1)

A.1

Here

u_{m}

(⁠

m = 1, \dots, M

⁠) is the concentration of molecular species

U_{m}

(⁠

m = 1, \dots, M

⁠) in system.

r_{j} (k_{j}; u)

is the reaction rate of reaction j

(j = 1, \dots, J)

which is depending on reaction parameter

k_{j}

and

u

⁠. From the reaction equation for reaction

j

⁠:

j : \sum_{m = 1}^{M} y_{m}^{j} U_{m} \to \sum_{m = 1}^{M} {\bar{y}}_{m}^{j} U_{m} (A .2)

A.2

ν_{m j}

is defined using the stoichiometric coefficient of the reaction as follows:

ν_{m j} := {\bar{y}}_{m}^{j} - y_{m}^{j} .

M \times J

matrix

ν = (ν_{m j})

is termed a stoichiometric matrix. A stoichiometric matrix is usually expressed using integer components.

In the following, we consider a case where the dimension of the left null space of

ν

(⁠

ker (ν^{T}) = {d \in R^{M} | d^{T} ν = 0}

⁠) is zero,

dim (ker (ν^{T})) = 0

⁠, i.e. there is no

d \in R^{M}

satisfying

d^{T} ν = 0

except for the zero vector. Here, using

rank (ν) = M - dim (ker (ν^{T})) = J - dim (ker (ν))

⁠,

ν

is full rank (⁠

rank (ν) = M \leq J

⁠). This corresponds to the case where there is no conserved quantity in the dynamical system. If

d

is not a zero vector, we multiply both sides of the system of differential equations by

d^{T}

from the left to obtain:

d^{T} \frac{d}{d t} u = d^{T} ν r = 0 (A .3)

A.3

from which we determine

d^{T} u = \sum_{m = 1}^{M} d_{m} u_{m}

is the quantity conserved in the dynamics.

Reaction rate $r_{j}$ is a function of the concentration of the substrate and the regulating molecule. However, the specific form of the rate function is unknown. Each reaction rate $r_{j}$ also has a parameter $k_{j}$ corresponding to the activity or amount of the respective enzyme.

Let us consider a steady state of this system (A.1). In steady state

(\hat{r}, \hat{u})

⁠, total inflow and total outflow is equivalent for each m, and we have:

\sum_{j = 1}^{J} ν_{m j} {\hat{r}}_{j} (k_{j}; \hat{u}) = 0, m = 1, \dots, M . (A .4)

A.4

Here, from the requirement that there exists a steady state reaction rate

r

that is not a zero vector,

dim (ker (ν)) > 0

holds for null space of

ν = (ν_{m j})

(right null space,

ker (ν) = {c \in R^{J} | ν c = 0}

⁠). The steady state reaction rate vectors

\hat{r} = ({\hat{r}}_{1}, \dots, {\hat{r}}_{J})^{T}

is expressed using:

\hat{r} = \sum_{n = 1}^{N} μ^{n} c^{n} (A .5)

A.5

Here,

c^{n} (n = 1, \dots, N)

is the basis vector of the null space of the stoichiometric matrix. Since the rank of stoichiometric matrix

ν

is M, the dimension of the right null space is

N = dim (ker (ν^{T})) = J - M

⁠. Basis vector

c^{n}

corresponds to the path of the steady flow determined from the network structure. The basis vector

c^{n}

can also be interpreted as the cycle created by the reaction pathway on the network, or the pathway connecting the inflow to the outflow of the system. The definition of a cycle here is different from that in graph theory: a cycle consists of only reactions (edges) and does not include vertices. It is also permissible to follow the reaction arrows in reverse, since the components of

c^{n}

can be negative. The path from the inflow to the outflow of a system can also be said to be a cycle in a broad sense by assuming a path through a hypothetical state outside the system (the inflow and outflow of the entire system correspond to the outflow and inflow from the outside state, respectively).

As can be seen from (A.4), the steady state solution cannot be determined unless

r_{j}

is given specifically. However, even without determining a specific steady state solution, if we assume the existence of a steady state solution, it is possible to determine how the steady state solution varies with respect to the change in enzyme activity

k_{j^{*}}

(k_{j^{*}} \to k_{j^{*}} + δ k_{j^{*}})

from the network structure alone [1–4]. Assume that the steady state holds both before and after change in enzyme activity

k_{j^{*}}

⁠. The change in reaction rate,

δ_{j^{*}} r_{j}

in response to change in activity

(k_{j^{*}} \to k_{j^{*}} + δ k_{j^{*}})

is given in the following form using the total derivative:

δ_{j^{*}} r_{j} = [\frac{\partial r_{j}}{\partial k_{j^{*}}} + \sum_{m = 1}^{M} \frac{\partial r_{j}}{\partial u_{m}} \frac{\partial u_{m}}{\partial k_{j^{*}}}] δ k_{j^{*}} (j, j^{*} = 1, \dots, J) . (A .6)

A.6

Note that reaction rate

r_{j} (k_{j}; u)

is a function of the reaction parameter

k_{j}

and molecular concentrations

u

⁠, and

u

may depend on

k_{j^{*}}

⁠. Below, all partial derivatives are evaluated at the steady state.

(A.6) implies that

δ_{j^{*}} r_{j}

is decomposed into a direct response to a change in

k_{j^{*}}

(first term), and an indirect response via a change in concentration

u_{m}

(second term).

\partial r_{j} / \partial k_{j^{*}}

is positive only when

j = j^{*}

⁠, and is zero when

j \neq j^{*}

⁠. On the other hand, from (A.5), the change in the reaction rate can be expressed as

δ_{j^{*}} r_{j} = \sum_{n = 1}^{N} δ_{j^{*}} μ^{n} c_{j}^{n}

⁠. Here,

δ_{j^{*}} μ^{n}

is the change in

μ^{n}

due to the change in enzymatic activity

(k_{j^{*}} \to k_{j^{*}} + δ k_{j^{*}})

⁠. Hence, we obtain:

\frac{\partial r_{j}}{\partial k_{j^{*}}} δ k_{j^{*}} + \sum_{m = 1}^{M} \frac{\partial r_{j}}{\partial u_{m}} \frac{\partial u_{m}}{\partial k_{j^{*}}} δ k_{j^{*}} = \sum_{n = 1}^{N} δ_{j^{*}} μ^{n} c_{j}^{n} (j, j^{*} = 1, \dots, J) . (A .7)

A.7

By a simple transformation, we have:

\sum_{m = 1}^{M} r_{j m} δ_{j^{*}} u_{m} - \sum_{n = 1}^{N} δ_{j^{*}} μ^{n} c_{j}^{n} = - \frac{\partial r_{j}}{\partial k_{j^{*}}} δ k_{j^{*}} (j, j^{*} = 1, \dots, J) . (A .8)

A.8

Here, we set

r_{j m} := \frac{\partial r_{j}}{\partial u_{m}}

and

δ_{j^{*}} u_{m} := \frac{\partial u_{m}}{\partial k_{j^{*}}} δ k_{j^{*}}

⁠. For an easier understanding, we show another form of Equation (A.8):

(\begin{matrix} r_{j 1} & \dots & r_{j m} & \dots & r_{j M} \end{matrix}) (\begin{matrix} δ_{j^{*}} u_{1} \\ ⋮ \\ δ_{j^{*}} u_{m} \\ ⋮ \\ δ_{j^{*}} u_{M} \end{matrix}) + (\begin{matrix} - c_{j}^{1} & \dots & - c_{j}^{n} & \dots & - c_{j}^{N} \end{matrix}) (\begin{matrix} δ_{j^{*}} μ^{1} \\ ⋮ \\ δ_{j^{*}} μ^{n} \\ ⋮ \\ δ_{j^{*}} μ^{N} \end{matrix}) = {\begin{matrix} - \frac{\partial r_{j}}{\partial k_{j^{*}}} δ k_{j^{*}} (j = j^{*}) \\ 0 (j \neq j^{*}) \end{matrix} (j, j^{*} = 1, \dots, J) . (A .9)

A.9

From (A.8) (or (A.9)) by taking j in row and

j^{*}

in column, we have:

(\begin{matrix} r_{11} & \dots & r_{1 M} \\ ⋮ & ⋮ \\ r_{J 1} & \dots & r_{J M} \end{matrix} | \begin{matrix} - c_{1}^{1} & \dots & - c_{1}^{N} \\ ⋮ & ⋮ \\ - c_{J}^{1} & \dots & - c_{J}^{N} \end{matrix}) (\begin{matrix} \begin{matrix} δ_{1} u_{1} & \dots & δ_{J} u_{1} \\ ⋮ & ⋮ \\ δ_{1} u_{M} & \dots & δ_{J} u_{M} \end{matrix} \\ \bar{\begin{matrix} δ_{1} μ^{1} & \dots & δ_{J} μ^{1} \\ ⋮ & ⋮ \\ δ_{1} μ^{N} & \dots & δ_{J} μ^{N} \end{matrix}} \end{matrix}) = - (\begin{matrix} \frac{\partial r_{1}}{\partial k_{1}} δ k_{1} & 0 \\ ⋱ \\ 0 & \frac{\partial r_{J}}{\partial k_{J}} δ k_{J} \end{matrix}) . (A .10)

A.10

Let

A := (r_{j m} | - c^{n})

be the first matrix on the left-hand side of Equation (A.10). The right-hand side is a diagonal matrix, where only the diagonal components

(j = j^{*})

take nonzero value. Since we are only interested in the qualitative response here, the right-hand side can be taken as an identity matrix. Assuming that A has an inverse matrix, we obtain:

S \equiv (\begin{matrix} \begin{matrix} δ_{1} u_{1} & \dots & δ_{J} u_{1} \\ ⋮ & ⋮ \\ δ_{1} u_{M} & \dots & δ_{J} u_{M} \end{matrix} \\ \bar{\begin{matrix} δ_{1} μ^{1} & \dots & δ_{J} μ^{1} \\ ⋮ & ⋮ \\ δ_{1} μ^{N} & \dots & δ_{J} μ^{N} \end{matrix}} \end{matrix}) \propto - A^{- 1} . (A .11)

A.11

In Equation (A.11), each column identifies perturbed reaction parameter

(j^{*} = 1, \dots, J)

and each row identify molecular species

(m = 1, \dots, M)

or basis vector

(n = 1, \dots, N)

of

ν

⁠. Each component of

S

shows the response of concentration

u_{m}

of each molecule m (or that of steady flow

μ^{n}

of cycle

n

⁠) to the change in enzyme activity

k_{j^{*}}

of reaction

j^{*}

⁠. Equation (A.11) gives the response to an increase in the activity (e.g. overexpression) of each enzyme. If one is interested in the response to a decrease in activity (e.g. knockdown), the minus sign on the right-hand side should be simply removed.

The above explanation shows the steps of derivation of the formula. To analyze the perturbation responses from the network structure in practice, one of the following criteria is used depending on the case.

(1) Determining the presence or absence of a response, i.e. zero or nonzero of $δ_{j^{*}} u_{m}$ or $δ_{j^{*}} μ^{n}$ ⁠.

Zero or non-zero of each component of the matrix

A

is determined solely from the structure of the reaction network, i.e. stoichiometric matrix and dependence of reaction rate functions. That is,

r_{j m}

is constitutively zero when reaction j does not depend on the molecular species m, and nonzero when it does (unless special values of concentrations are chosen). The components of

c^{n}

are also determined from the stoichiometric matrix. In other words, the distribution of the zero (nonzero) components of matrix

A

reflects the structure of the reaction network. The presence or absence of change in each concentration or each reaction rate to perturbation of each reaction parameter, that is, zero or nonzero of each component in the matrix

S

is determined only from the distribution of the zero component (non-zero component) in the matrix

A

⁠. That are determined only from the structure of the reaction network.

The Equation (A.7) uses the partial derivative, it appears to hold only when the change $δ k_{j^{*}}$ in enzymatic activity $k_{j^{*}}$ is small. However, the constitutive zero components of matrix $A$ is determined from the network structure alone, and does not depend on constant $k$ or the concentration $u$ ⁠. Thus, the pattern of responses calculated by $A^{- 1}$ (the distribution of zero components in matrix $S$ ⁠) will be always the same, even after repeated application of the change in $k_{j^{*}}$ ⁠, $(k_{j^{*}} \to k_{j^{*}} + δ k_{j^{*}})$ ⁠. In other words, the result does not change, even if the change $δ k_{j^{*}}$ in parameter is large.

(2) Determining the direction of the response, i.e. sign of $δ_{j^{*}} u_{m}$ or $δ_{j^{*}} μ^{n}$ ⁠.

Suppose we want to determine not only the presence or absence of a response, but also direction of responses (whether the concentration or reaction rate increases or decreases). In this case, we need to assume that the reaction rates

r_{j}

are monotonic functions of the concentration of the influencing molecular species. For example, assume that it is an increasing function of substrate molecules and an increasing (positive regulation) or decreasing (negative regulation) function of regulator molecules. As a result,

r_{j m}

is determined to be zero, positive, or negative. The matrix

A

with sign information in addition to zero or nonzero information of each component is given. Then we calculate inverse of the matrix

A

to determine the sign of

δ_{j^{*}} u_{m}

and

δ_{j^{*}} μ^{n}

⁠. Empirically, unless the structure of the network is too complex, it is possible to determine the sign of the response uniquely in many cases. When networks are complicated, though zero or nonzero of components are always determined uniquely, it may happen that the signs of nonzero components are not uniquely decided. This is because the components in the inverse matrix contain multiple

r_{j m}

's and their signs may be different. In this case, the signs of the responses can be determined again by giving the relative magnitudes of some

r_{j m}

's contained in the components.

As mentioned at the beginning, we have assumed that the dimension of the left null space of

ν

is 0, that is, there is no conserved quantity in the system. If there are conserved quantities in the system, the basis vectors

d^{l}

(l = 1, \dots, L)

of the left null space

ker (ν^{T}) = {d \in R^{M} | d^{T} ν = 0}

(here

L = dim (ker (ν^{T}))

of stoichiometric matrix

ν

can be used to construct the extended matrix

A

⁠:

A := (\begin{matrix} \underline{\begin{matrix} r_{11} & \dots & r_{1 M} \\ ⋮ & ⋮ \\ r_{J 1} & \dots & r_{J M} \end{matrix} | \begin{matrix} - c^{1} & \dots & - c^{N} \end{matrix}} \\ \begin{matrix} - {(d^{1})}^{T} \\ ⋮ \\ - {(d^{L})}^{T} \end{matrix} | \begin{matrix} 0 & \dots & 0 \\ ⋮ & ⋮ \\ 0 & \dots & 0 \end{matrix} \end{matrix}) . (A .12)

A.12

By calculating the inverse of this matrix

A

⁠, the response to perturbation is determined:

S \equiv (\begin{matrix} \underline{\begin{matrix} δ_{1} u_{1} & \dots & δ_{J} u_{1} \\ ⋮ & ⋮ \\ δ_{1} u_{M} & \dots & δ_{J} u_{M} \end{matrix} | \begin{matrix} δ_{J + 1} u_{1} & \dots & δ_{J + L} u_{1} \\ ⋮ & ⋮ \\ δ_{J + 1} u_{M} & \dots & δ_{J + L} u_{M} \end{matrix}} \\ \begin{matrix} δ_{1} μ^{1} & \dots & δ_{J} μ^{1} \\ ⋮ & ⋮ \\ δ_{1} μ^{N} & \dots & δ_{J} μ^{N} \end{matrix} | \begin{matrix} δ_{J + 1} μ^{1} & \dots & δ_{J + L} μ^{1} \\ ⋮ & ⋮ \\ δ_{J + 1} μ^{N} & \dots & δ_{J + L} μ^{N} \end{matrix} \end{matrix}) \propto - A^{- 1} . (A .13)

A.13

Here,

δ_{J + l} u_{m}

is the change in the steady state concentration of molecule m when the conserved quantity

\sum_{m^{'} = 1}^{M} d_{m^{'}}^{l} u_{m^{'}} = T^{l}

(determined from the initial value) corresponding to the

l

th basis vector

d^{l}

of the left null space of

ν

is perturbed, and

δ_{J + l} μ^{n}

is the change in the steady state flow along the basis vector

c^{n}

in the null space of

ν

when

T^{l}

is perturbed. For details of the analysis, see [4].

Appendix B: Law of localization

For the chemical reaction system and the A matrix discussed in the Appendix A, the Law of Localization and the buffering structure are given as follows [3,4].

Definition (Buffering structure)

Consider an arbitrary subgraph

Γ = (M^{*}, J^{*})

on a chemical reaction network. Here,

M^{*}

is a subset of molecules and

J^{*}

is a subset of reactions. For the null space of the stoichiometric matrix

v

⁠, let

N^{*} (J^{*})

be the set of linearly independent basis vectors such that the components are zero except for the reactions contained in

J^{*}

⁠. Here we introduce the following index:

χ (Γ) \equiv | M^{*} | - | J^{*} | + | N^{*} (J^{*}) | . (B .1)

B.1

Subgraph

Γ

is called a ‘buffering structure' when

Γ

satisfies the following conditions:

(i) (All reactions dependent on molecules in $M^{*}$ ⁠) $\subseteq J^{*}$
(ii) $χ (Γ) = 0$

Note 1: Condition (i) implies that, when considering subgraphs, ‘Once a subset

M^{*}

of molecules is determined, all reactions that use them as substrates or regulators are included in the subset

J^{*}

of reactions’. After satisfying this condition, arbitrary reactions may be added to

J^{*}

⁠.

Note 2: $| N^{*} (J^{*}) |$ is the number of linearly independent basis vectors contained in $J^{*}$ within null space $ν$ ⁠, i.e. the number of independent cycles consisting only of reactions contained in $J^{*}$ ⁠. Thus, $χ (Γ)$ indicates ‘(number of molecular species) − (number of reactions) + (number of cycles).' Since it corresponds to the alternating sum of the number of zero-, one-, and two-dimensional structures in the subgraph $Γ$ ⁠, it is strongly related to the Euler characteristic. Note that the value of $χ (Γ)$ is generally negative or zero when the matrix $A$ is regular.

Theorem (Law of Localization)

Consider a steady state of a chemical reaction system. When a reaction parameter in a buffering structure is changed, the steady-state concentration of molecules outside the buffering structure and the steady-state rate of reactions outside the buffering structure do not change. In other words, the response of the system remains inside the buffering structure.

Proof

When a buffering structure exists in a chemical reaction system, by exchanging rows and columns of matrix

A

⁠, we can set a block upper triangular matrix such that:

A = (\begin{matrix} A_{Γ} & * \\ 0 & A_{\bar{Γ}} \end{matrix}) . (B .2)

B.2

Here,

A_{Γ}

is a square submatrix of size

| J^{*} | \times (| M^{*} | + | N^{*} (J^{*}) |)

⁠, where

| J^{*} | = | M^{*} | + | N^{*} (J^{*}) |

⁠. Each row of

A_{Γ}

corresponds to a reaction in

J^{*}

⁠, each column corresponds to a molecular species in

M^{*}

⁠, or to a basis vector of

ker (ν)

in

N^{*} (J^{*})

⁠.

A_{Γ}

and

A_{\bar{Γ}}

are not zero matrices, and when

A

is regular, they are also invertible.

0

is a submatrix of size

(| J | - | J^{*} |) \times (| M^{*} | + | N^{*} (J^{*}) |)

such that all components are zero. From the distribution of the zero components in the matrix, the inverse matrix has the same distribution of zero components, i.e.

S \propto - A^{- 1} = (\begin{matrix} A_{Γ}^{- 1} & * \\ 0 & A_{\bar{Γ}}^{- 1} \end{matrix}) . (B .3)

B.3

Again,

0

is a submatrix of size

(| J | - | M^{*} | - | N^{*} (J^{*}) |) \times | J^{*} |

such that the components are all zero. Each component of

0

indicates the change in concentration of an external molecule (⁠

m \notin M^{*}

⁠), or the change in an external steady-state flow

n \notin N^{*} (J^{*})

when a reaction parameter inside the buffering structure

j^{*} \in J^{*}

is perturbed. Thus, when a chemical reaction system contains a buffering structure, the responses to a perturbation of a reaction parameter inside the buffering structure remain inside the buffering structure. (End of proof)

Competing Interests

The authors declare that there are no competing interests associated with the manuscript.

References

1

Mochizuki

,

A.

and

Fiedler

,

B.

(

2015

)

Sensitivity of chemical reaction networks: a structural approach. 1. Examples and the carbon metabolic network

.

J. Theor. Biol.

367

,

189

–

202

https://doi.org/10.1016/j.jtbi.2014.10.025

Google Scholar

Crossref

PubMed

2

Fiedler

,

B.

and

Mochizuki

,

A.

(

2015

)

Sensitivity of chemical reaction networks: a structural approach. 2. Regular monomolecular systems

.

Math. Meth. Appl. Sci.

38

,

3381

–

3600

https://doi.org/10.1002/mma.3436

Google Scholar

Crossref

3

Okada

,

T.

and

Mochizuki

,

A.

(

2016

)

Law of localization in chemical reaction networks

.

Phys. Rev. Lett.

117

,

048101

https://doi.org/10.1103/PhysRevLett.117.048101

Google Scholar

Crossref

PubMed

4

Okada

,

T.

and

Mochizuki

,

A.

(

2017

)

Sensitivity and network topology in chemical reaction systems

.

Phys. Rev. E

96

,

022322

https://doi.org/10.1103/PhysRevE.96.022322

Google Scholar

Crossref

PubMed

5

Ishii

,

N.

,

Nakahigashi

,

K.

,

Baba

,

T.

,

Robert

,

M.

,

Soga

,

T.

,

Kanai

,

A.

et al (

2007

)

Multiple high-throughput analyses monitor the response of E. coli to perturbations

.

Science

316

,

593

–

597

https://doi.org/10.1126/science.1132067

Google Scholar

Crossref

PubMed

6

Crown

,

S.B.

and

Antoniewicz

,

M.R.

(

2013

)

Publishing ¹³C metabolic flux analysis studies: a review and future perspectives

.

Metab. Eng.

20

,

42

–

48

https://doi.org/10.1016/j.ymben.2013.08.005

Google Scholar

Crossref

PubMed

7

Clarke

,

B.L.

(

1975

)

Theorems on chemical network stability

.

J. Chem. Phys.

62

,

773

https://doi.org/10.1063/1.430524

Google Scholar

Crossref

8

Clarke

,

B.L.

(

1975

)

Stability of topologically similar chemical networks

.

J. Chem. Phys.

62

,

3726

https://doi.org/10.1063/1.430969

Google Scholar

Crossref

9

Chellaboina

,

V.

,

Bhat

,

S.P.

,

Haddad

,

W.M.

and

Bernstein

,

D.S.

(

2009

)

Modeling and analysis of mass-action kinetics

.

IEEE Cont. Sys. Mag.

29

,

60

–

78

https://doi.org/10.1109/MCS.2009.932926

Google Scholar

Crossref

10

Feinberg

,

M.

(

1995

)

The existence and uniqueness of steady states for a class of chemical reaction networks

.

Arch. Rational Mech. Anal.

132

,

311

–

370

https://doi.org/10.1007/BF00375614

Google Scholar

Crossref

11

Craciun

,

G.

and

Feinberg

,

M.

(

2006

)

Multiple equilibria in complex chemical reaction networks: I. The injectivity property

.

SIAM J. App. Math.

65

,

1526

–

1546

https://doi.org/10.1137/S0036139904440278

Google Scholar

Crossref

12

Craciun

,

G.

and

Feinberg

,

M.

(

2006

)

Multiple equilibria in complex chemical reaction networks: II. The species-reactions graph

.

SIAM J. App. Math.

66

,

1321

–

1338

https://doi.org/10.1137/050634177

Google Scholar

Crossref

13

Shinar

,

G.

and

Feinberg

,

M.

(

2013

)

Concordant chemical reaction networks and the species-reaction graph

.

Math. Biosci.

241

,

1

–

23

https://doi.org/10.1016/j.mbs.2012.08.002

Google Scholar

Crossref

PubMed

14

Kacser

,

H.

and

Burns

,

J.A.

(

1973

)

The control of flux

.

Symp. Soc. Exp. Biol.

27

,

65

–

104

.

Reprinted in Biochem. Soc. Trans. 23, 341–366, 1995

https://doi.org/10.1042/bst0230341

Google Scholar

Crossref

PubMed

15

Heinrich

,

R.

and

Rapoport

,

T.A.

(

1974

)

A linear steady-state treatment of enzymatic chains; general properties, control and effector strength

.

Eur. J. Biochem.

42

,

89

–

95

https://doi.org/10.1111/j.1432-1033.1974.tb03318.x

Google Scholar

Crossref

PubMed

16

Fell

,

D.A.

(

1992

)

Metabolic control analysis: a survey of its theoretical and experimental development

.

Biochem. J.

286

,

313

–

330

https://doi.org/10.1042/bj2860313

Google Scholar

Crossref

PubMed

17

Stephanopoulos

,

G.N.

,

Aristidou

,

A.A.

and

Nielsen

,

J.

(

1998

)

Metabolic Engineering: Principles and Methodologies

,

Academic Press

,

San Diego, California, USA

Google Scholar

18

Ferjani

,

A.

,

Kawade

,

K.

,

Asaoka

,

M.

,

Oikawa

,

A.

,

Okada

,

T.

,

Mochizuki

,

A.

et al (

2018

)

Pyrophosphate inhibits gluconeogenesis by restricting UDP-glucose formation in vivo

.

Sci. Rep.

8

,

14696

https://doi.org/10.1038/s41598-018-32894-1

Google Scholar

Crossref

PubMed

19

Steuer

,

R.

,

Waldherr

,

S.

,

Sourjik

,

V.

and

Kollmann

,

M.

(

2011

)

Robust signal processing in living cells

.

PLoS Comput. Biol.

7

,

e1002218

https://doi.org/10.1371/journal.pcbi.1002218

Google Scholar

Crossref

PubMed

20

Hens

,

C.

,

Harush

,

U.

,

Haber

,

S.

,

Cohen

,

R.

and

Barzel

,

B.

(

2019

)

Spatiotemporal signal propagation in complex networks

.

Nat. Phys.

15

,

403

–

412

https://doi.org/10.1038/s41567-018-0409-0

Google Scholar

Crossref

21

Okada

,

T.

,

Tsai

,

J.-C.

and

Mochizuki

,

A.

(

2018

)

Structural bifurcation analysis in chemical reaction networks

.

Phys. Rev. E

98

,

012417

https://doi.org/10.1103/PhysRevE.98.012417

Google Scholar

Crossref

PubMed

22

Okada

,

T.

,

Mochizuki

,

A.

,

Furuta

,

M.

and

Tsai

,

J.-C.

(

2021

)

Flux-augmented bifurcation analysis in chemical reaction network systems

.

Phys. Rev. E

103

,

062212

https://doi.org/10.1103/PhysRevE.103.062212

Google Scholar

Crossref

PubMed

2022

This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY-NC-ND).

A structural approach to understanding enzymatic regulation of chemical reaction networks

Introduction

Chemical reaction networks

Reaction network of the central metabolic pathway of bacteria (Escherichia coli).

Dynamics of chemical reaction networks

Sensitivity analysis of a small reaction network.

Structural sensitivity analysis

Example of chemical reaction system.

Application (1): plant metabolism

Structural sensitivity analysis of the plant metabolism system.

Law of localization

Examples of buffering structures.

Application (2): central metabolic pathway

The perturbation response network of the E. coli central metabolic pathway obtained from the analysis.

Discussion

Appendix A: Mathematics of structural sensitivity analysis

Appendix B: Law of localization

Competing Interests

References

Cited By

Get Email Alerts

CONNECT

EXPLORE

Cover Image

A structural approach to understanding enzymatic regulation of chemical reaction networks

Introduction

Chemical reaction networks

Reaction network of the central metabolic pathway of bacteria (Escherichia coli).

Dynamics of chemical reaction networks

Sensitivity analysis of a small reaction network.

Structural sensitivity analysis

Example of chemical reaction system.

Application (1): plant metabolism

Structural sensitivity analysis of the plant metabolism system.

Law of localization

Examples of buffering structures.

Application (2): central metabolic pathway

The perturbation response network of the E. coli central metabolic pathway obtained from the analysis.

Discussion

Appendix A: Mathematics of structural sensitivity analysis

Appendix B: Law of localization

Competing Interests

References

Cited By

Get Email Alerts

CONNECT

EXPLORE

This Feature Is Available To Subscribers Only