Theoretical approaches for the dynamics of complex biological systems from information of networks

Abstract

Modern biology has provided many examples of large networks describing the interactions between multiple species of bio-molecules. It is believed that the dynamics of molecular activities based on such networks are the origin of biological functions. On the other hand, we have a limited understanding for dynamics of molecular activity based on networks. To overcome this problem, we have developed two structural theories, by which the important aspects of the dynamical properties of the system are determined only from information on the network structure, without assuming other quantitative details. The first theory, named Linkage Logic, determines a subset of molecules in regulatory networks, by which any long-term dynamical behavior of the whole system can be identified/controlled. The second theory, named Structural Sensitivity Analysis, determines the sensitivity responses of the steady state of chemical reaction networks to perturbations of the reaction rate. The first and second theories investigate the dynamical properties of regulatory and reaction networks, respectively. The first theory targets the attractors of the regulatory network systems, whereas the second theory applies only to the steady states of the reaction network systems, but predicts their detailed behavior. To demonstrate the utility of our methods several biological network systems, and show they are practically useful to analyze behaviors of biological systems.

1. Introduction

Modern biology has provided us with numerous examples of large networks describing the interactions between multiple species of molecules, such as genes and proteins (e.g.^1)–4)). The dynamics of such network-based molecular activities are widely considered as the origin of biological functions. For example, the circadian rhythms observed in many species are produced by periodic oscillations of gene activities.⁵⁾ Diversities of differentiated cells are considered to be caused by the diversity of steady states of gene expressions.⁶⁾ A major objective in modern biology is to understand biological functions in terms of the dynamics of the activity of bio-molecules, based on experimentally determined information of networks.

Biological networks, which graphically represent the interactions between bio-molecules, are divided into several classes. A “regulatory network” describes the influences between bio-molecules, i.e., each arrow in the graph indicates “who controls who”. A “reaction network” describes the state transitions undergone by the reactions, i.e., each arrow in the graph indicates “who becomes who”. Examples of regulatory networks and reaction networks are abundant in databases. They have been derived from the accumulated experimental results in molecular genetics and biochemistry. On the other hand, the dynamics resulting from, and encoded in, such complex network systems are not sufficiently understood.

One remaining difficulty is the observation of dynamic processes. The activity dynamics of bio-molecules are not easily observed with sufficient time resolution. Most of the data obtained by present experimental methods are snapshots of molecular activities rather than time tracks. The second problem is the reliability of the network information. At present, the networks compiled from many studies of biological systems may be incomplete, because of the complexity and working cost of experimentally identifying the reaction edges. The third and largest problem is that the dynamics cannot be determined from the network information alone. The regulatory edges provide only qualitative information on the dependencies among activities of bio-molecules in the system. They lack essential quantitative details, such as the regulatory/reaction functions, parameter values of reaction rates, and initial states.

To overcome these problems, we have developed two mathematical theories that derive the dynamical properties of complex biological systems from information of the regulatory/reaction linkages alone.^7)–11) The first theory, called Linkage Logic, is developed for “regulatory networks”. Linkage Logic applies a broad class of ordinary differential equation (ODE) systems. The only mathematical requirements are dissipativity (bounded solutions) and a decay condition (See [2]). It ensures that:

1. i)    All non-transient dynamical behaviors of the whole system can be identified faithfully by the measurement of only a subset of variables in the system, and
2. ii)    The subset is determined from the regulatory linkage alone as a feedback vertex set (FVS) of the network graph, as explained below.

The second theory is developed for “reaction networks”. Some mathematical approaches study dynamical properties of reaction systems from the network information.^12),13) One possible experimental method to reaction networks examines their sensitivity: each enzyme mediating a reaction in the system is increased/decreased or knocked out separately, and the responses in concentrations of chemicals or their fluxes are measured by mass spectrometry. There are some mathematical studies to analyze sensitivity of chemical reaction networks. Kacser and Burns¹⁴⁾ and Heinrich and Rapoport¹⁵⁾ 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. However, without sufficient analyses, these studies did not find any response patterns related to topology of networks, nor any general laws connecting responses and topology. We will present a different, and simpler, mathematical framework which enables us to calculate sensitivity of large systems qualitatively. We analyze sensitivity by a structural method and determine relations between topology of networks and sensitivity responses.

In our second theory, called structural sensitivity analysis, we study the concentration (and flux) changes in a reaction system under perturbation of reaction rate parameters, assuming that the system is in a steady state. The theory determines the qualitative change in the chemical concentrations and fluxes in the system solely from information on the reaction network. We did not assume any specific functions modelling the reactions, or their reaction parameters. Initially, we assume the monotonicity of reaction function of substrates (see [4]). We then relax this assumption and use information of dependence of reaction function only (see [4′]).

In this paper, we demonstrate the power of our theories in understanding the dynamics of practical biological systems from experimentally determined networks. In particular, we discuss and check the consistency between the network structures and the observed dynamics of molecular activities, while lacking quantitative information on the unobserved or hidden parts of the dynamics.

2. Linkage Logic

Let us explain Linkage Logic, our first structural theory of regulator networks. The molecular activity dynamics are modeled as the following system of ODEs.   

\begin{equation} \frac{dx_{k}}{dt} = f_{k}(\boldsymbol{{x}}_{I_{k}}) - d_{k}(x_{k})\quad k = 1, \cdots,N \end{equation} [1]

Here, x_k denotes the activity of bio-molecule k at time t, and f_k is a positive non-linear function describing the activity enhancement of molecule k (which may be called a regulatory function), and d_k is a positive increasing function describing the activity decay.^18),19) The set I_k ⊆ {1, ⋯, N} is a subset of molecules that regulate molecule k. In other words, I_k is the input set of k. The boldface notation $\boldsymbol{{x}}_{I_{k}}$ denotes the vector of components x_j with j ∈ I_k. The set I_k includes k, i.e., k ∈ I_k, if and only if the self-activation of molecule k exceeds its self-repression and decay.

More generally, we broaden the classes of the ODE system of regulatory networks:   

\begin{equation} \frac{dx_{k}}{dt} = F_{k}(x_{k},\boldsymbol{{x}}_{I_{k}})\quad k = 1, \cdots,N \end{equation} [2]

Here $F_{k}(x_{k},\boldsymbol{{x}}_{I_{k}})$ are any non-linear functions. The dynamics are assumed as dissipative i.e., every solution is bounded after a positive time. The final assumption is the decay condition $\partial _{1}F_{k}(x_{k},\boldsymbol{{x}}_{I_{k}}) < 0$, where $\partial_{1} $ denotes the first partial derivative with respect to the first argument. This system allows for possible self-activation. If it exists, the self-activation of molecule k is expressed as k ∈ I_k. Note that Eq. [2] includes decay-rate regulation by molecules other than k.^9),10) The regulatory network Γ simply consists of vertices (all biomolecular “species” k), and directed edges (all j → k such that j is an element of the input set I_k of k). In the following, we consider the dynamics of system [1] or [2] based on the information of the regulatory network I_k, only, without using any further information.^8)–10)

Here we explain and combine two concepts from different mathematical fields. The first concept is the feedback vertex set (FVS) from graph theory. A FVS is a subset of vertices in a directed graph whose removal leaves the graph without directed cycles.²⁰⁾ To clarify the concept, examples of small networks with the highlight of feedback vertex sets are illustrated in Fig. 1.

Fig. 1.

Examples of small regulatory networks. The directed edges show the regulatory interactions between the nodes. The gray vertices highlight a selected minimal feedback vertex set in each case (a)–(e). Modified from Mochizuki et al. (2013).¹⁰⁾

The second concept is the determining node from dynamical systems theory.^21)–23) In the setting of [1] and [2], we call a subset of variables J ⊆ {1, ⋯, N} determining nodes, if and only if the convergence of variables in J, for any two trajectories, implies the convergence of all variables of these trajectories. More precisely, let x(t) and $\tilde{\boldsymbol{{x}}}(t)$ be two solutions of [1] or [2]. Then $\boldsymbol{{x}}(t) - \tilde{\boldsymbol{{x}}}(t) \to 0$ holds in the large-time limit t → +∞, for all components k ∈ {1, ⋯, N}, if the two solutions satisfy $\boldsymbol{{x}}_{J}(t) - \tilde{\boldsymbol{{x}}}_{J}(t) \to 0$, for all components k in the subset J ⊆ {1, ⋯, N}. In other words, if the long-term dynamics of the determining nodes J are given, then the dynamics of the whole system are uniquely determined in the long-term.

We mathematically proved that (i) any feedback vertex set of a regulatory network is a set of determining nodes of the dynamics on the network. Conversely, (ii) if a vertex set is determining, for all choices of nonlinearities compatible with the network structure, then it is a feedback vertex set.

The first statement (i) is proved in straightforward manner. We can easily show that $\boldsymbol{{x}}(t) - \tilde{\boldsymbol{{x}}}(t) \to 0$ holds for t → +∞, if two solutions satisfy $\boldsymbol{{x}}_{J}(t) - \tilde{\boldsymbol{{x}}}_{J}(t) \to 0$ on a feedback vertex set J of a regulatory network. The second statement (ii) is proved by showing contrapositive.⁹⁾ Here, we provide a brief intuitive explanation of our theory. We first consider a single regulation in a network. Of course, if the dynamics of the input vertices are given, the long-term dynamics of the downward vertex are uniquely determined. If the regulatory function leading to the lower vertex, the dynamics is not determined constructively, but is still determined uniquely. Next, let us consider a system of a regulatory network including several vertices and edges as in Fig. 2. Inductively repeating our previous argument downward through the network, we can uniquely determine the whole system dynamics if the dynamics of an appropriate subset of vertices is given. Of course, the dynamics of the total system can be uniquely determined only when that subset is appropriately chosen; all of the remaining vertices should be downward from the vertices in that subset. Third, let us consider a problem: how can we minimize the subset on which the dynamics are given? The solution is the minimal feedback vertex set I. Indeed, if I is a feedback vertex set of a graph with vertex set Γ = {1, 2, ⋯, N}, then by definition of the feedback vertex set I, all remaining vertices in K = Γ ∖ I can be ordered downward from I. This implies that the dynamics of all vertices in the system can be uniquely determined by induction from the dynamics on the feedback vertex set I.

Fig. 2.

Intuitive interpretation of the theory. If the dynamics are known at the gray vertices, the dynamics of the remaining vertices are uniquely determined. The set of vertices on which the dynamics are given can be reduced to a minimal feedback vertex set (in this case, the single vertex marked by a red circle). Details are given in the text. Modified from Mochizuki et al. (2013).¹⁰⁾

3. Cell differentiation in Ascidiacea development

We consider the gene regulatory network determined by Imai et al. (2006) (Fig. 3a),³⁾ which governs cell-differentiation in the development of the marine filter feeder Ciona intestinalis (Ascidiacea) from the 16-cell stage to the tail-bud stage. In the focal period, the gene activities between the cells diversify over time. The tail-bud stage is characterized by 13 different gene expression patterns, depending on the cells’ position in the body. The tail-bud stage, with its regional gene expression patters, continues for a relatively longer period than the previous developmental stages. The system is expected to be sufficiently flexible to produce many steady states of gene activities corresponding to the differentiated cell types.

Fig. 3.

Gene regulatory network of Ascidiacea development. (a) Network based on Imai et al. (2006).³⁾ The original network includes 16 genes with self-repression. We removed these repressive self-loops, because self-repression is subsumed under degradation in our formulation [2]. (b) Reduced network after successive removal of nodes without input or output. The reduced network contains a minimal feedback vertex set with a single vertex, FoxD-a/b. Modified from Mochizuki et al. (2013).¹⁰⁾

As a preparation, we removed the vertices that do not receive any regulations or do not regulate any vertices. Each of these top genes converges to a stationary point and the bottom genes provide the system outputs; therefore, they make no contribution to the diversity of attractors. After removing the irrelevant vertices, the network reduces to the 7-vertex network of Fig. 3b. The network possesses a minimal feedback vertex set comprising a single feedback vertex, FoxD-a/b. All long-term dynamics on the global attractor possibly generated by this gene regulatory network can therefore be identified by measurement of the activity of the single gene FoxD-a/b, if the network information is correct. In present developmental biology studies, gene expressions are usually interpreted in a discrete and binary manner, i.e., are classified as active (1) or inactive (0). Of course, it is impossible to identify the diversity of the 13 different gene expressions by only two points of one binary variable.

The result can be interpreted in two ways. The first possibility is that the diversity of differentiated cells may be captured not by binary, but by continuous values of activity of FoxD-a/b. In this case, if experimental biologists could measure FoxD-a/b activity by a more precise method, the diversity of cell states might be indicated by the different activity levels of FoxD-a/b. The second interpretation is that the present understanding of the gene regulatory network of Ascidiacea is not sufficient to explain the diversity of gene expressions at the tail bud stage. In other words, the observed cell differentiation in Ascidiacea might be largely governed by unknown regulatory linkages.

4. Structural sensitivity analysis

We now introduce our second structural theory, developed for reaction networks. Dubbed structural sensitivity analysis, this theory structurally determines the sensitivity responses of a chemical reaction network to perturbations of any reaction in the system. We study the change in concentrations in a reaction system under perturbation of reaction rate parameter, assuming that the system is in a steady state (Mochizuki & Fiedler, 2015).¹¹⁾ In the following, chemicals and reactions are labeled as m (m = 1, …, M) and i (i = 1, …, R), respectively. In general, a system state is specified by the participant concentrations x_m(t), which obey the following differential equations   

\begin{equation} \frac{dx_{m}}{dt} = \sum_{i = 1}^{R}\nu_{mi}W_{i}(k_{i};x) \end{equation} [3]

The matrix ν in [3] is called a stoichiometric matrix, and the function W_i is called a flux, which depends on the metabolite concentrations and the reaction rate parameter k_i. The reaction rate specifies the amount/activity of enzyme mediating the reaction. The flux functions are assumed to increase with substrate concentration, but are not constrained to a specific form. Thus we set   

\begin{align} &\frac{\partial W_{i}}{\partial x_{m}} > 0\ \text{if $x_{m}$ is a substrate of reaction $i$; }\\ &\text{otherwise}\ \frac{\partial W_{i}}{\partial x_{m}} = 0. \end{align} [4]

Hereafter, $\frac{\partial W_{i}}{\partial x_{m}}$ is abbreviated as r_im.

In this framework, an enzyme knockdown of the j-th reaction corresponds to a change in the corresponding reaction rate parameter; that is, k_j → k_j + δk_j. We assume equilibrium of this system both before knockdown and after knockdown leading the following condition:   

\begin{align} \sum_{i = 1}^{R}\nu_{mi}\delta_{j}W_{i} &= \sum_{i = 1}^{R}\nu_{mi}\left(\frac{\partial W_{i}}{\partial k_{j}} + \sum_{m' = 1}^{M}\frac{\partial W_{i}}{\partial x_{m'}}\frac{dx_{m'}}{dk_{j}} \right)\delta k_{j} \\ &= 0 \end{align} [5]

for any m. Here, $\delta _{j}\vec{W}$ is the flux change induced by the parameter change.

As shown in Mochizuki and Fiedler (2015),¹¹⁾ the changes in steady state concentration $\delta _{j}\vec{x} \equiv \frac{d\vec{x}}{dk_{j}}\delta k_{j}$ and flux $\delta _{j}\vec{W}$ under the perturbation k_j → k_j + δk_j are solely determined from the network structure. The following matrix computation simultaneously obtains the chemicals’ responses to perturbations of each reaction rate:   

\begin{equation} \left(\frac{\delta_{1}\vec{x} \cdots \delta_{R}\vec{x}}{\delta_{1}\vec{\mu} \cdots \delta_{R}\vec{\mu}} \right) \propto (\mathbf{A}^{-1})_{mj} \equiv \mathbf{S} \end{equation} [6]

where the square matrix A is given as   

\begin{equation} \mathbf{A} = (r_{im}| - \vec{c}_{1} \cdots - \vec{c}_{N_{k}}) \end{equation} [7]

and the inverse of A is called the sensitivity matrix S. The flux change is then determined as follows:   

\begin{equation} (\delta_{1}\vec{W} \cdots \delta_{R}\vec{W}) = (\vec{c}_{1} \cdots \vec{c}_{N_{k}})(\delta_{1}\vec{\mu} \cdots \delta_{R}\vec{\mu}) \end{equation} [8]

where $\vec{c}_{n}$ is the basis of the kernel space of the stoichiometric matrix ν, and δ_jμⁿ are its coefficients. Note that the distribution of nonzero entries in A-matrix reflects the structure of the reaction network. Therefore, we can determine the qualitative response and flux of each chemical from the distribution of nonzero entries in A-matrix only. This implies that our theory depends only on the structure of the reaction network. To integrate regulatory effects such as allosteric control, we can generalize our method to [4′] by relaxing [4]:   

\begin{align} &\frac{\partial W_{i}}{\partial x_{m}} \neq 0\ \text{if $x_{m}$ influences reaction $i$, }\\ &\text{otherwise}\ \frac{\partial W_{i}}{\partial x_{m}} = 0. \end{align} [4′]

This modification places additional nonzero entries in A-matrix. The sensitivity is again determined by [6].

5. Examples

Let us demonstrate our structural sensitivity analysis using simple examples.

5-1. Single pathway.

First, we consider the single reaction pathway shown in Fig. 4. A-matrix and the sensitivity matrix S are given by   

\begin{align} \mathbf{A} & = \begin{pmatrix} 0 & 0 & 0 & -1\\ r_{2A} & 0 & 0 & -1\\ 0 & r_{3B} & 0 & -1\\ 0 & 0 & r_{4C} & -1 \end{pmatrix} \\ \mathbf{S} & = \begin{pmatrix} -r_{2A}^{-1} & r_{2A}^{-1} & 0 & 0\\ -r_{3B}^{-1} & 0 & r_{3B}^{-1} & 0\\ -r_{4C}^{-1} & 0 & 0 & r_{4C}^{-1}\\ -1 & 0 & 0 & 0 \end{pmatrix} \end{align} [9]

The result shows that the flux changes if, and only if, the top reaction 1 (input to the system) is perturbed. Then all chemical concentrations in the system change, accordingly. However, perturbations of reactions 2, 3 or 4 alter only the input reactant chemical of the perturbed reaction; the other chemical concentrations and fluxes are unchanged. This result is easily interpreted as follows: as the reaction rate of a chemical increases, the concentration of the upstream chemical reactant decreases to compensate the change in the reaction rate. This buffering effect prevents the perturbation from propagating upward or downward from its original locus. The equilibria of the single pathway were calculated in numerical simulations of several reaction functions. Our mathematical result was consistently matched by the simulation results.

Fig. 4.

Sensitivity analysis of a single reaction pathway of chemical reactions with a feed 1 and an exit reaction 4. (a–d) Changes in concentrations and fluxes induced by perturbations of reaction rates k₁ to k₄ from top (a) to bottom (d). Red triangles indicate the perturbed reactions j = 1, ⋯, 4. The plus or minus signs adjacent to the circles indicate an increase or decrease in the concentrations of chemicals A, B, C, respectively. Red bold circles (and arrows) indicate increases in concentrations (and fluxes). Red dashed circles (and arrows) indicate decrease in concentrations (and fluxes). Modified from Mochizuki and Fiedler (2015).¹¹⁾

5-2. Feedback circuit.

As a second example, we analyze the small network shown in Fig. 5. The network consists of 6 reactions and 4 chemicals, including one feedback loop via one side branch. The matrix A-matrix is given as:   

\begin{equation} \mathbf{A} = \begin{pmatrix} 0 0 0 0 -1 0\\ r_{2A} 0 0 0 -1 -1\\ 0 r_{3B} 0 0 -1 -1\\ 0 0 r_{4C} 0 0 -1\\ 0 0 0 r_{5D} 0 -1\\ 0 0 r_{6C} 0 -1 0 \end{pmatrix} \end{equation} [10]

The inverse of the A-matrix gives sensitivity, responses of chemical concentrations and of fluxes to the perturbations as:   

\begin{equation} {\mathbf{S} = \begin{pmatrix} -(r_{4C} + r_{6C})/r_{2A}r_{6C} r_{2A}^{-1} 0 -r_{2A}^{-1} 0 r_{4C}/r_{2A}r_{6C}\\ -(r_{4C} + r_{6C})/r_{3B}r_{6C} 0 r_{3B}^{-1} -r_{3B}^{-1} 0 r_{4C}/r_{3B}r_{6C}\\ -r_{6C}^{-1} 0 0 0 0 r_{6C}^{-1}\\ -r_{4C}/r_{5D}r_{6C} 0 0 -r_{5D}^{-1} r_{5D}^{-1} r_{4C}/r_{5D}r_{6C}\\ -1 0 0 0 0 0\\ -r_{4C}/r_{6C} 0 0 -1 0 r_{4C}/r_{6C} \end{pmatrix}} \end{equation} [11]

In [11], the columns indicate the perturbed reactions j = 1, ⋯, 6, and the rows indicate the sensitivity responses of the chemical concentrations A, B, C, D (first four rows) and fluxes (last two rows). Again, only a change in the top feed reaction w₁ can simultaneously affect all chemicals and fluxes (see Fig. 5a). Reactions j = 2, 3 and 5 are buffered by their input reactants A, B, D, respectively, and their perturbations propagate no further through the network (see Fig. 5b, c, e).

Fig. 5.

Sensitivity analysis of a simple network with one feedback loop. (a–f) Perturbed reactions (indicated by red triangles) are j = 1, ⋯, 6. Changes in concentrations and fluxes induced by perturbations of reaction rates k₁, …, k₆ are indicated. See also the legend of Fig. 1. Modified from Mochizuki and Fiedler (2015).¹¹⁾

Let us compare the results of perturbations to reaction j = 4 (fourth column, Fig. 5d) and j = 6 (last column, Fig. 5f). A decrease in the reaction constant of reaction 4 affects a large part of the network (Fig. 5d). We understand the result: the chemical concentrations D, A, B along the reaction cycle 4, 5, 2, 3 downward of the perturbed reaction j = 4 simply decrease, along with the fluxes of these reactions. Perturbations of reaction j = 6 also largely disturb the network (Fig. 5f). This time, however, the changes in concentrations do not appear in the chemical concentrations further downstream, but in the side branch cycle 4, 5, 2, 3 of the perturbed reaction. Also, the signs of the concentration and flux changes are reversed in Fig. 5d and 5f. In other words, the perturbations into reaction branches j = 4 and 6, both of which symmetrically emanate from chemical C, yield asymmetric sensitivity responses.

As shown in Fig. 5f, a perturbation to reaction j = 6 increases the concentration of chemical C. This is a compensatory response that maintains the output flow (w₆ of reaction 6) at the level of the unchanged feed flow w₁ (restriction along the vertical pathways). On the other hand, any change in chemical C spreads to the second branch, reaction 4, which emanates from C. The consequent increase in x_C increases the flux and concentrations in the feedback loop, thus buffering the output flow w₆ against the externally applied perturbation.

6. Metabolic networks

As a real life application, we analyze the chemical reaction system of carbon metabolism in Escherichia coli. For this purpose, we adopt the network presented in Ishii (2007)⁴⁾ with minor modifications (specifically, we consider 28 metabolites and 43 reactions). The reaction network is illustrated in Fig. 6. We constructed the A-matrix of [7] and inverted it to obtain the sensitivity matrix S of [6].

Fig. 6.

Reaction network of the TCA cycle in the carbon metabolic network of E. coli. Colors summarize the influence patterns of the flux responses shown in Fig. 7, with some details omitted. The flux influence patterns are summarized as: i) yellow→yellow, ii) red→{red, yellow}, iii) brown→{red, yellow}, iv) light blue→light blue, v) blue→{blue (including dashed blue), light blue}, vi) green→{green (including dashed), blue, light blue, brown, red, yellow}, vii) black→{black (including dashed), green, blue, light blue, brown, red, yellow}. Modified from Mochizuki and Fiedler (2015).¹¹⁾

The colors of the fluxes in Fig. 7 indicate the hierarchy in the response patterns of the fluxes. The nonzero responses of flux to perturbations are simply summarized as follows. The fluxes are categorized into 6 color classes linked by transitive hierarchical relations; namely, black→green, green→{red, yellow, blue}, red→yellow, blue→light blue. The global carbon metabolism of E. coli is insensitive to the reaction sets of the pentose phosphate pathway (red, yellow) and the TCA cycle (blue, light blue). Perturbations to a reaction in these sets influence all members in the own sets, but exert no influence on other sets. On the other hand, perturbing some reactions in the glycolysis pathway (green) influences almost all reactions in the network. We emphasize that these insights are solely gained by our mathematical analysis of the sensitivity matrix S. Other than the network structure, no a priori biological knowledge of the metabolism is required to reach these conclusions.

Fig. 7.

Hierarchy of nonzero response patterns of the carbon metabolism network shown in Fig. 6. Modified from Mochizuki and Fiedler (2015).¹¹⁾

7. Discussion

We introduced two structural theories that determine the important aspects of the dynamical behaviors of complex systems from the information of network structures alone. We also applied the theories to real biological systems. Modern biology has acquired a rich body of information describing the interactions among bio-molecules in networks. However, to specify the network dynamics, we require qualitative details that are currently difficult to obtain. We believe that the proposed theories will play an important role in analyzing the behaviors of complex biological systems in the near future.

In complex regulatory networks, we identified small FVSs whose measurements would assuredly identify the dynamical behavior of the whole system. Dynamic behaviors include steady states, periodic oscillations and quasi-periodic oscillations. In some examples, the networks may be insufficient to explain the diversity of the observed dynamics in actual biological systems. Our theory also provides a rational criterion for selecting the key molecules to control the system; that is, by prescribing the FVS dynamics, we can sufficiently control the whole system.

For reaction networks, we developed our structural sensitivity analysis, which determines the rate sensitivity of the systems in a structural manner without assuming a functional form. In this analysis, the qualitative responses are solely determined from the structure of the network. Specific choices of the reaction functions and their rate parameters are not required. The characteristic responses of the system depend on the network structure and the location of the rate perturbation in the network. We also analyzed the patterns of nonzero flux responses, and discovered a clear hierarchy in the flux response.

Our function-free structural approaches make predictions from the network structures. On the other hand, at present, our knowledge of regulatory/reaction networks in many organisms is probably incomplete. Combined with specific experimental measurements, our theory might provide a useful tool for revealing unknown reactions or regulations of network systems. For example, we may calculate the sensitivity of hypothetical modifications that introduce or omit some reactions or regulations, and compare the calculations with the results of specific experiments. If modifying the “known” network considerably enhances agreement with experimental results, then the modified network would likely represent the real system. Worsening disagreement, however, would immediately falsify the network, annulling the possibility of more suitable reaction functions or fudgy rate constants.

As demonstrated in this paper, our theories can elucidate the dynamics of complex biological systems. Combined with actual data of molecular activities, our theory could become a powerful tool in molecular biology for predicting the existence of unknown molecules or regulatory linkages. In addition, we hope that our theory will support the rational selective acquisition of currently missing biological data, and further elucidate the complex and fascinating mechanisms of life.

Profile

Atsushi Mochizuki was born in Shizuoka, Japan, in 1969. He graduated from the Faculty of Sciences, Kyoto University, in 1994, and obtained his Ph.D. in 1999 from Kyushu University. He was promoted to assistant professor in 1998 at Kyushu University, and to associate professor in 2002 at National Institute for Basic Biology. Since then, he has been director of his own research group. He was promoted to chief scientist in 2008 at RIKEN. He has studied biological phenomena using mathematical and computational methods to understand the principles of complex dynamics in biology. His research achievements include theoretical studies for regulatory networks, sensitivity of chemical reaction systems, pattern formations in development, evolution of gene regulation, and so on. His final goal is to establish a new bioscience that progresses by the repeats of theoretical predictions and experimental verifications near future. For his achievements, he received 11th JSPS Prize in 2014.

Acknowledgements

This work was supported in part by the JST PRESTO program, CREST program of Japan, an HFSP program grant. We express our sincere thanks to the following people for their useful comments: T. Akutsu, G. Craciun, B. Fiedler, N. Go, Y. Iwasa, S. Ishihara, K. Ishimatsu, D. Kiga, H. Kishino, H. Kokubu, H. Kurata, S. Mano, H. Matano, Y. Nishiura, H. Oka, K. Okada, T. Okada, G. Rangel, N. Shigesada, H. Stuke and M. Wakayama.

References

• 1)   Davidson,  E.H.,  Rast,  J.P.,  Oliveri,  P.,  Ransick,  A.,  Calestani,  C.,  Yuh,  C.H.,  Minokawa,  T.,  Amore,  G.,  Hinman,  V.,  Arenas-Mena,  C.,  Otim,  O.,  Brown,  C.T.,  Livi,  C.B.,  Lee,  P.Y.,  Revilla,  R.,  Rust,  A.G.,  Pan,  Z.J.,  Schilstra,  M.J.,  Clarke,  P.J.C.,  Arnone,  M.I.,  Rowen,  L.,  Cameron,  R.A.,  McClay,  D.R.,  Hood,  L. and  Bolouri,  H. (2002) A genomic regulatory network for development. Science 295, 1669–1678.
• 2)   Oda,  K.,  Matsuoka,  Y.,  Funahashi,  A. and  Kitano,  H. (2005) A comprehensive pathway map of epidermal growth factor receptor signaling. Mol. Syst. Biol. 1, 0010.
• 3)   Imai,  K.S.,  Levine,  M.,  Satoh,  N. and  Satou,  Y. (2006) Regulatory blueprint for a chordate embryo. Science 312, 1183–1187.
• 4)   Ishii,  N.,  Nakahigashi,  K.,  Baba,  T.,  Robert,  M.,  Soga,  T.,  Kanai,  A.,  Hirasawa,  T.,  Naba,  M.,  Hirai,  K.,  Hoque,  A.,  Ho,  P.Y.,  Kakazu,  Y.,  Sugawara,  K.,  Igarashi,  S.,  Harada,  S.,  Masuda,  T.,  Sugiyama,  N.,  Togashi,  T.,  Hasegawa,  M.,  Takai,  Y.,  Yugi,  K.,  Arakawa,  K.,  Iwata,  N.,  Toya,  Y.,  Nakayama,  Y.,  Nishioka,  T.,  Shimizu,  K.,  Mori,  H. and  Tomita,  M. (2007) Multiple high-throughput analyses monitor the response of E. coli to perturbations. Science 316, 593–597.
• 5)   Mirsky,  H.P.,  Liu,  A.C.,  Welsh,  D.K.,  Kay,  S.A. and  Doyle,  F.J. III (2009) A model of the cell-autonomous mammalian circadian clock. Proc. Natl. Acad. Sci. U.S.A. 106, 11107–11112.
• 6)   Albert,  R. and  Othmer,  H.G. (2003) The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J. Theor. Biol. 223, 1–18.
• 7)   Mochizuki,  A. (2008) Structure of regulatory networks and diversity of gene expression patterns. J. Theor. Biol. 250, 307–321.
• 8)   Mochizuki,  A. and  Saito,  D. (2010) Analyzing steady states of dynamics of bio-molecules from the structure of regulatory networks. J. Theor. Biol. 266, 323–335.
• 9)   Fiedler,  B.,  Mochizuki,  A.,  Kurosawa,  G. and  Saito,  D. (2013) Dynamics and control at feedback vertex sets I: Informative and determining nodes in regulatory networks. J. Dyn. Differential Eqns. 25, 563–604.
• 10)   Mochizuki,  A.,  Fiedler,  B.,  Kurosawa,  G. and  Saito,  D. (2013) Dynamics and control at feedback vertex sets. II: A faithful monitor to determine the diversity of molecular activities in regulatory networks. J. Theor. Biol. 335, 130–146.
• 11)   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.
• 12)   Feinberg,  M. (1995) The existence and uniqueness of steady states for a class of chemical reaction networks. Arch. Ration. Mech. Anal. 132, 311–370.
• 13)   Craciun,  G. and  Feinberg,  M. (2006) Multiple equilibria in complex chemical reaction networks: The species-reactions graph. SIAM J. Appl. Math. 66 (4), 1321–1338.
• 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.
• 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.
• 16)   Fell,  D.A. (1992) Metabolic control analysis: a survey of its theoretical and experimental development. Biochem. J. 286, 313–330.
• 17)  Stephanopoulos, G.N., Aristidou, A.A. and Nielsen, J. (1998) Metabolic Engineering: Principles and methodologies. Academic Press, San Diego.
• 18)   Glass,  L. and  Kauffman,  S.A. (1973) The logical analysis of continuous, non-linear biochemical control networks. J. Theor. Biol. 39, 103–129.
• 19)   Mochizuki,  A. (2005) An analytical study of the number of steady states in gene regulatory networks. J. Theor. Biol. 236, 291–310.
• 20)   Akutsu,  T.,  Kuhara,  S.,  Maruyama,  O. and  Miyano,  S. (1998) A system for identifying genetic networks from gene expression patterns produced by gene disruptions and overexpressions. Genome Informatics. 9, 151–160.
• 21)   Foias,  C. and  Temam,  R. (1984) Determination of the solutions of the Navier-Stokes equations by a set of nodal values. Math. Comput. 43, 117–133.
• 22)   Foias,  C. and  Titi,  E.S. (1991) Determining nodes, finite differences schemes and inertial manifolds. Nonlinearity 4, 135–153.
• 23)  Takens, F. (2010) Reconstruction theory and nonlinear time series analysis. In Handbook of Dynamical Systems 3 (eds. Broer, H. et al.). Elsevier, Amsterdam, pp. 345–377.