Roles of network topology in the relaxation dynamics of simple chemical reaction network models

In the following sections, we study the link between the network structure and the relaxation behavior of the accompanying ODE. First, we generate ca. 2000 minimum reaction networks (see Methods for the parameters we used) and check if there is a steady-state attractor. Almost all the instances of the reaction system has deficiency more than one34,35,36,37,38. However, as far as we have confirmed, all the networks each have a unique steady-state attractor. We have confirmed that the unistability holds even if the reaction rate constants are randomly chosen, the Michaelis-Menten type rate equation is utilized, and non-minimal networks are used for constructing kinetic models (see Fig. S1-4). This robustness of unistability might be related to the absence of feedback regulations in the system such as substrate-level regulation and catalytic reactions.For probing the relaxation characteristics, we applied perturbations to the attractor. Since the concentrations of the chemicals span a wide range at the attractors, we apply perturbations in a multiplicative manner. In addition, to normalize the strength of the perturbations, we generated the initial points on the surface of the \(|{\mathscr {C}}_0|\)-dimensional hypersphere, centred at the steady-state attractor, with a given diameter D. In concrete terms, a single initial point \({\varvec{x}}_\textrm{ini}\) is given by \(\ln {\varvec{x}}_\textrm{ini}=\ln {\varvec{x}}_\textrm{att}+D\cdot {\varvec{r}}/\Vert {\varvec{r}}\Vert\) with \({\varvec{x}}_\textrm{att}\) is the steady-state attractor of the model, and \({\varvec{r}}\) is a vector of the random numbers each lying in the interval [0, 1], respectively (see Materials and Methods section for more details).For a sufficiently small value of D, such that the dynamics could be described adequately by linearizing the ODE around the steady-state, it is known that the relaxation dynamics would be exponential, with a rate determined by the eignevalues of the Jacobi matrix. Hence, we choose a value of D that is large enough to go beyond this linear regime so we can explore how the non-linearities of the ODEs result in non-exponential relaxation.We simulated the accompanying ODE (Eq. (6)) of ca. 2000 minimum networks with \(N_\textrm{ptb}(=256)\) initial points for each network. Visually, we found three typical types of relaxation: exponential, plateau, and power-law type. Examples of each are shown in Fig.1C–F. There appeared to be three groups in terms of classification of the networks: (1) the networks exhibiting only exponential relaxation regardless of the initial concentrations, (2) the networks showing only the plateau type relaxation regardless of the initial concentrations, and (3) the networks displaying either the plateau or the power-law type of relaxation depending on the initial concentrations. At this stage, we do not know if the plateau exhibited by the second and the third type of networks are different in any sense, but for later arguments we refer to them as “Metastable Plateau” (independent of initial conditions) and “Confined Plateau” (dependent on initial conditions) in Fig. 1D and E. “Confined” here means two things: the first meaning is that the initial conditions leading to the plateau behavior are confined in a certain region in the phase space. The second meaning is the confinement of the trajectory in the so-called “stoichiometric cone”. See the subsection III for more details.To make the classification more quantitative, we computed two quantities: the relaxation time and the migration length. Here we define the relaxation time as the time point at which the Euclidean distance, in logarithmic scale, between the state and the steady-state attractor first becomes smaller than a certain threshold value \(\epsilon\), that is$$\begin{aligned} T_\textrm{relax}=\min \{t\ |\ \Vert \ln {\varvec{x}}(t)-\ln {\varvec{x}}_{att}\Vert <\epsilon \}. \end{aligned}$$
(7)
The choice of this definition is for the nature of the relaxations in the non-linear regime. In the linear response regime, the relaxation time is characterized by the rate constant of the exponential relaxation, while as will be shown, the model shows power-law relaxation where no such trivial time-scale exists. As discussed in39, the choice of \(\epsilon\) is crucial in this type of definition. While we adopted \(\epsilon =10^{-6}\), it is checked that our conclusions are unchanged as long as \(\epsilon\) value is in the range \((10^{-6},10^{-2})\).Next, to distinguish between the plateau and the power-law relaxation types, we quantified the total length of the trajectory in concentration space during a given time period. It is measured by the “migration length” defined as$$\begin{aligned} L= \int _{a}^{b}\Bigl \Vert \frac{{\textrm{d}}\ln {\varvec{x}}}{dt}\Bigr \Vert dt, \end{aligned}$$
(8)
where we choose \(a=100\) and \(b=0.01/\phi\) in order to quantify the “motion” of the trajectory in the intermediate time scale (sufficiently larger than the time scale of the reaction rate constants, which are unity, and sufficiently smaller than that of the degradation, \(1/\phi\)). The migration length is computed in the logarithmic scale to take the effect of low-concentration chemicals into account. Additionally, we compute \(L_\textrm{max}\) given by Eq. (8) with \(b = 10^{10}\) (the end of the simulation) as the upper bound of the integral.Our prior expectation, from visual examination of the trajectories (e.g., Fig.1C–F), was that since plateau type trajectories are almost static in the time range \(t\in [100,0.01/\phi ]\), the normalized migration length \(L/L_\textrm{max}\) would have a relatively lower value, compared to that for power-law type dynamics. We also expected that exponential relaxation dynamics will result in lower \(T_\textrm{relax}\) values than the other two types.Examples of the distribution of migration length of two sets of networks are presented in Fig. 1G. The upper panel is the distribution over all initial conditions for one network that shows Metastable-plateau dynamics, and the bottom panel is for one network exhibiting both confined-plateau and power-law relaxation. The bottom distribution shows a clear double-peak due to bimodality of the relaxation dynamics, while the upper distribution is rather unimodal. We suppose that if networks exhibit a bimodal distribution due to the presence of both power-law relaxation in addition to confined-plateau relaxation, the average value of the normalized migration length \(L/L_\textrm{max}\) (averaged over all initial conditions) would be larger than the corresponding value for metastable-plateau networks. Thus, we utilize \(\langle L/L_\textrm{max}\rangle\) as an indicator of the variety of the relaxation behavior.The average relaxation time and the average migration length are computed for each network, and we plotted the distribution over the networks. As shown in Fig. 1H, the distribution has four distinct peaks, and we classified the networks from type-I to type-IV according to the nearest peak. While the separation of the type-II and type-IV peaks is less clear, the distribution-based classification of the networks matched very well to the visual grouping of the models: the type-I networks (small \(T_\textrm{relax}\) and large \(L/L_\textrm{max}\)) exhibit only exponential relaxation; most of the type-II networks showed metastable-plateau relaxation; and, the type-III networks and most type-IV networks exhibited both the confined-plateau and the power-law type relaxation depending on the initial concentrations.Note that the final steady-state attractor is reached at the latest by \(t\sim 1/\phi\), which is the timescale of the spontaneous degradation, and the qualitative differences between the different types of relaxation only appear in the \(t\ll 1/\phi\) region. These observations mean that the difference in the relaxation behavior during \(t\ll 1/\phi\) highlight structural features of the reaction networks and would not be changed by varying the value of \(1/\phi\) as long as it is sufficiently large. What structural features of the networks result in the four different relaxation types observed in Fig.1H? We now computationally show that the relaxation behavior is indeed determined by three features of the stoichiometric matrix (which in turn encodes the structure of the network), namely, (I) the rank gap, (II) dimension of the left null-space, and (III) the stoichiometric cone.Fig. 1(A) The fully-connected network for \(M=2\) and \(L=6\). (B) Two examples of the minimum network which synthesizes a single target from two sources. The sources and target chemicals are indicated by the red arrows and the blue arrow, respectively. In (A and B), the bipartite representation of the chemical reaction network is adopted. For example, the reaction \(S_1+S_2\rightarrow P\) is depicted as follows (see inset of the panel B); first, the chemical nodes representing \(S_1\) and \(S_2\) (the gray disks) are wired to the reaction node (the pink diamond), and then, an edge connects the reaction node and the chemical node of the chemical P (the gray disk). Directions of the reactions are chosen to be consistent with the steady-state flux distribution. (C–F) Four typical relaxation dynamics emerged from the minimum networks: Exponential, Metastable-Plateau, Confined-Plateau, and Power-law relaxation. The vertical and horizontal axes are on a logarithmic scale. The four dynamics are classified into three category based on the networks types that exhibit the corresponding dynamics. (G) Example distributions of the migration length (see text for definition). The top panel is the distribution over all initial conditions for one network showing metastable-plateau type relaxation, and the bottom panel is from a network showing both the plateau and power-law type relaxations. (H) Classification of the minimum networks based on the average relaxation time (see text for definition) and the average migration length.The rank gapThe distinct relaxation behaviors are observed over timescales shorter than \(1/\phi\), and thus, should stem from the dynamics (Eq. (6)) without the degradation term. The clear difference between the exponential relaxation and the others is the relaxation time. The exponential type trajectories already reach the attractor well before t exceeds \(1/\phi\). We now show how measuring the ranks of the stoichiometric matrix and some enlargements can inform us about the reachability to the attractor in the absence of the degradation term.First, we introduce additional notation. Let \(S_1\) denote the enlarged matrix \((S_0|{\varvec{U}}^{(0)}|{\varvec{U}}^{(1)}|\cdots |{\varvec{U}}^{(N_n-1)})\), where \({\varvec{U}}^{(j)}\)’s are vectors of the stoichiometry of the source chemical uptake reactions, i.e., \(U_i^{(j)}\) is unity if the ith chemical is the jth source chemical and zero otherwise. By vertically stacking \({\varvec{v}}_0\) and \({\varvec{u}}\) in Eq. (6) and denoting it as \({\varvec{v}}_1\), we can rewrite Eq. (6) with \(S_1\) instead of \(S_0\). Then, removing the degradation term, we obtain$$\begin{aligned} \dot{\varvec{x}}=S_1{\varvec{v}}_1-{\varvec{p}}. \end{aligned}$$
(9)
Note that Eq. (9) represents the dynamics of chemical concentration driven by the source chemical uptake, internal reactions, and the target synthesis. This equation approximates the original equation Eq. (6) for \(t\ll 1/\phi.\)For there to exist a steady-state of this equation, there must be a pair of vectors, \({\varvec{v}}_1^*\) and \({\varvec{p}}^*\), satisfying \(S_1{\varvec{v}}_1^*={\varvec{p}}^*\). Since we are not interested in the trivial solution \({\varvec{v}}_1^*={\varvec{0}}\), the necessary and sufficient condition for the existence of a steady-state is that the rows of \(S_1\) and \({\varvec{P}}\) (The stoichiometric part of the harvest reaction) are linearly dependent, i.e, \(\textrm{rank} S_1=\textrm{rank} (S_1|{\varvec{P}})\), where \({\varvec{P}}\) is the stoichiometry of the target synthesis (i.e., the ith element is unity if the ith chemical is one of the target precursors and zero for the rest).We found that the rank gap \(\delta \equiv \textrm{rank} (S_1|{\varvec{P}})-\textrm{rank}S_1\) for all type-I networks (showing only the exponential relaxation) is zero. Interestingly, all type-II networks also satisfied the \(\delta =0\) condition. On the other hand, all the type-III and type-IV networks have \(\delta =1\), i.e., in the absence of degradation no steady-state exists.This implies that there are two types of plateau relaxation. The plateaux exhibited by type-II networks are the steady-states of the model equation without degradation (Eq. (6)). However, the plateaux exhibited by type-III and IV networks are not the steady-state of Eq. (6), but are extremely slow transient dynamics towards the final steady-state for the model with degradation. Indeed, we find that the type-II plateaux correspond to what we had earlier termed “metastable-plateau”, while the type-III and IV plateux correspond to what we called “confined-plateau”.The left null-spaceThe rank gap does not distinguish between type-I and type-II networks (which both have \(\delta =0\)). Fortunately, these networks exhibit a unique type of relaxation regardless of the initial concentrations (exponential or metastable-plateau). Thus, we can expect that the difference stems from the network structure.Figure 2 shows two time courses of the same type-II network but starting from different initial concentrations. The plateau regions differ not only in the concentrations of the chemicals but also in the rank order of the chemical concentrations (the identical color is used for each chemical species in both figure panels).Fig. 2Two time courses emerged from the same minimum type-II network when started with different initial concentrations. Both the vertical and horizontal axes are plotted in the logarithmic scale.Considering that each plateau corresponds to the steady-state of \(\dot{\varvec{x}}=S_1{\varvec{v}}_1-{\varvec{p}}\) and the model eventually reaches a unique final attractor, possible mechanisms to generate the variations are limited. We hypothesized that such type-II networks might be have combinations of chemical concentrations that, over timescales much shorter than \(1/\phi\), are constrained to be almost-conserved. If different initial conditions have different values for such conserved quantities, the system might evolve towards a different steady-state of \(\dot{\varvec{x}}=S_1{\varvec{v}}_1-{\varvec{p}}\) for each of these. We therefore computed the dimension of the left null-space of the enlarged stoichiometry matrix, \(I:=\textrm{dim}\left( \textrm{coker}\right) S_2\), where \(S_2=(S_1|{\varvec{P}})\), because the left null-space of the reaction stoichiometry represents the conserved quantities of the reaction dynamics. We found that the left null-space of all the type-I models was \(\{0\}\). In contrast, all the type-II models had non-zero I. Further, we found that the type-III models had no conserved quantity (\(I=0\)), while the type-IV networks had at least one (\(I\ge 1\)). (In the present minimum reaction models, there is a lower limit of I for the models with \(\textrm{rank}S_1=\textrm{rank}S_2\), given by \(M-N_n\), where M is the number of monomer types and \(N_n\) is the number of source chemicals . For more detail, see SI Text).So far, we have fully characterized the type-I and type-II networks. The type-I networks have no rank gap, \(\delta =\textrm{rank} S_1-\textrm{rank} S_2=0\), and no conserved quantity \(I=\textrm{dim}\left( \textrm{coker}\ S_2\right) =0\). Type-II networks also have no rank gap but have at least one conserved quantity, \(I\ge 1\).The stoichiometric coneNow we move onto the networks with \(\delta =1\) (Since the target formation reaction is represented as the single reaction (Eq. (5)), the rank gap between \(S_1\) and \(S_2\) is one at most). Due to the non-zero rank gap, there is no steady-state until t becomes comparable with \(1/\phi\), i.e., the spontaneous degradation starts affecting the dynamics of the concentration. Thus, in these type-III and type-IV networks, the relaxation dynamics over timescales smaller than \(1/\phi\) is not to some meta-stable almost-steady-state, but a slow transient dynamics towards the final attractor with degradation. With power-law relaxation there is a slow but continuous approach towards the final attractor, while with the confined-plateau dynamics the trajectory is practically static for a long time before suddenly “jumping” to the final attractor.We therefore hypothesized that in the power-law case, the attractor is reachable from the initial point, but the degradation term \(\phi {\varvec{x}}\) is necessary for making the final concentrations a stable steady-state. By contrast, in the confined-plateau case getting close to the attractor may itself be hindered in the absence of the degradation term. To check this, we introduce a subset of the phase space defined as$$\begin{aligned} SC({\varvec{x}}_\textrm{ini})=\{{\varvec{x}}\in {\mathbb R}^{|{\mathscr {C}}_0|}_+|{\varvec{x}}={\varvec{x}}_\textrm{ini}+S_2{\varvec{v}}_2,{\varvec{v}}_2\in {\mathbb R}^{|{\mathscr {R}}_0|},J^{(\textrm{tgt})}\ge 0\}. \end{aligned}$$\(SC({\varvec{x}}_\textrm{ini})\) is the region in the phase space “reachable” from the initial point by controlling the chemical reaction fluxes in any way except reversely harvesting the target chemicals. Note that in the definition, no flux balance (steady-state condition) is required. Due to the non-negativity condition (\(J^{({\varvec{t}gt})}\ge 0\)) and the non-zero rank gap \(\delta\), \(SC({\varvec{x}}_\textrm{ini})\) does not coincide with the whole space \({\mathbb R}^{|{\mathscr {C}}_0|}_+\). This \(SC({\varvec{x}}_\textrm{ini})\) is equivalent to the stoichiometric cone in the field of chemical reaction network theory5, and thus, we term \(SC({\varvec{x}})\) the stoichiometric cone corresponding to the initial condition \({\varvec{x}}\) in the following arguments.We can now mathematically formulate our hypothesis on the relationship between reachability and relaxation behaviors: for the initial points whose stoichiometric cone \(SC({\varvec{x}}_\textrm{ini})\) contains the final attractor, relaxation becomes power-law-like, whereas, if the stoichiometric cone of the initial point does not include the attractor, plateaux appear during the relaxation (for the algorithm of judging the inclusion relationship, see Methods).To examine the relationship between the relaxation type and the position of the initial concentrations in the phase space, we computed the average of the normalized migration length \(L/L_\textrm{max}\) separately for the initial concentrations whose stoichiometric cone does and does not contain the final attractor (Since the type-IV networks have conserved quantities, the stoichiometric cone of an initial point never contains the attractor unless the values of the conserved quantities at the initial point and the attractor are equal. This situation never happens by applying random perturbations to generate the initial points. Thus, we first search for the nearest point from the initial point, \(\tilde{{\varvec{x}}}_\textrm{ini}\), that has the same values as the conserved quantities with the attractor (see Materials and Methods). Then, we checked if the attractor is in the stoichiometric cone of \(\tilde{{\varvec{x}}}_\textrm{ini}\)). Figure 3 is the scatter plot of the two average values for all the minimum networks of type-III and IV that we constructed. The migration length averaged over the initial concentrations with \({\varvec{x}}_\textrm{att}\in SC({\varvec{x}}_\textrm{ini})\), \(\langle L/L_\textrm{max}\rangle _\textrm{in}\), and with \({\varvec{x}}_\textrm{att}\notin SC({\varvec{x}}_\textrm{ini})\), \(\langle L/L_\textrm{max}\rangle _\textrm{out}\) are plotted in the horizontal and vertical axis, respectively. In all networks, \(\langle L/L_\textrm{max}\rangle _\textrm{in}\) is larger than \(\langle L/L_\textrm{max}\rangle _\textrm{out}\). This result supports our intuition that some trajectories are confined in a region of the phase space and cannot get closer to the attractor.The difference between the initial conditions leading to the power-law and plateau relaxation is intuitively illustrated by a simple toy model shown in Fig. 4. The model equation is given by$$\begin{aligned} \frac{{\textrm{d}}{\varvec{x}}}{{\textrm{d}}t} = \begin{pmatrix} 1 & -2 \\ 1 & -1 \\ \end{pmatrix} \begin{pmatrix} C-x^{(A)}x^{(B)} \\ (x^{(A)})^2x^{(B)} \\ \end{pmatrix} -\phi {\varvec{x}}, \end{aligned}$$
(10)
where we set the direction of the reaction from left to right on Fig. 4, and C is the concentration of the source chemical (the leftmost chemical in the schematic). In the model, the second reaction is the target harvesting reaction.The first reaction produces one molecule each, but the second reaction consumes two molecules of chemical A and one molecule of chemical B. As a consequence of the imbalance, the stoichiometric cone of the model is given by$$\begin{aligned} SC({\varvec{x}}_\textrm{ini})=\{{\varvec{x}}\in {\mathbb R}^2_+|x^{(B)}-x^{(A)}\ge x_\textrm{ini}^{(B)}-x_\textrm{ini}^{(A)}\}. \end{aligned}$$This means that the quantity \(x^{(B)}-x^{(A)}\) only increase by the two reactions. In the model, we need to consume two molecules of A to consume one molecule of B. However, to increase the concentration of A, one needs to utilize the first reaction (note that the second reaction is irreversible), increasing the concentration of B by the same amount simultaneously.
Therefore, if the initial gap of the concentration is larger than the gap at the attractor, \(x_\textrm{att}^{(B)}-x_\textrm{att}^{(A)}<x_\textrm{ini}^{(B)}-x_\textrm{ini}^{(A)}\), the system cannot get closer than a certain distance for \(t\ll 1/\phi\). For instance, if \(x_\textrm{ini}^{(B)}-x_\textrm{att}^{(B)}>x_\textrm{ini}^{(A)}\) holds, the closest reacheable concentration of the chemical B is \(x_\textrm{ini}^{(B)}-x_\textrm{ini}^{(A)} > x_\textrm{att}^{(B)}\). In the contrast, if the stoichiometric cone contains the attractor, \(x_\textrm{att}^{(A)}\ll x_\textrm{ini}^{(A)}\) and \(x_\textrm{att}^{(B)}\gg x_\textrm{ini}^{(B)}\) for example, the system can approach arbitrary close to the attractor even in the absence of the degradation term.Typically, the quiescent dynamics on the confined plateau occur because the concentrations of one or a few chemicals become pretty low, as the chemical A in the toy example. Therefore, the concentrations at the plateau are not necessarily the closest possible state to the attractor in the stoichiometric cone but close to one of the axes of the Cartesian coordinate.Fig. 3The normalized migration length is averaged separately for the time courses with the initial points whose SC, defined in the text, contains (vertical axis) and does not contain (horizontal axis) the attractor. The \(y=x\) line is plotted as a guide for the eye. The distributions of \(\langle L/L_\textrm{max}\rangle _\textrm{in}\) and \(\langle L/L_\textrm{max}\rangle _\textrm{out}\) are shown along the horizontal and the vertical axis, respectively.Fig. 4A toy model to illustrate the power-law and plateau relaxation. The leftmost chemical is the external chemical whose concentration is fixed to a constant value.Fig. 5The summary of the relaxation type and network-structural features. Network types are classified according to the rank gap \(\delta\) and the dimension of the cokernel of the stoichiometry matrix \(S_2\), I. In the networks with \(\delta =1\), confined plateau or power-law dynamics emerges depending on whether the final attractor is in the stoichiometric cone of the initial point.Summary of the relaxation-structure relationshipTo summarize the above classification: Networks that have zero rank gap and no conserved quantity (type-I) always exhibit exponential relaxation dynamics. Networks that have zero rank gap but some conserved quantities (type-II) always exhibit plateau relaxation. Type III and type IV networks both have a rank gap of 1, but the former have no conserved quantities, while the latter have at least one. Both these networks exhibit either power-law relaxation or plateaux, depending on the initial conditions. For both types of networks, the type of relaxation for a given initial condition is determined by whether the final attractor is reachable from that initial point. If it is, then we observe power-law relaxation, and if it is not, then we observe confined plateau relaxation. Thus, from the stoichiometric matrices associated with a given network, by computing three network-structural quantities, namely, (1) the rank gap \(\delta\), (2) the dimension of the left null space \(I=\textrm{dim}\left( \textrm{coker}\right) S_2\), and (3) the stoichiometric cone, we can predict whether the relaxation will be exponential, power-law or plateau-type. A graphical summary of this result is presented in Fig. 5.