### Mathematical solution of the model

#### Degree distribution

First, we show an analytical solution for the degree distribution of the model via mean-field based analysis [38–40]. This analysis is based on a mean-field approximation, in which the many-body problem is considered as the one-body problem, and is widely used in the area of statistical mechanics of complex networks. Using the mean-field analysis, we can easily get the analytical solutions.

We here consider the time evolution of *k*_{
i
}, which is the degree (the number of edges) of node *i*. The degree of node *i* increases by one with the probability 1/*N*, where *N* is the total number of nodes, when Event I (a new metabolite and new reaction) occurs. When Event II occurs, two existing nodes are selected, and their degrees increase respectively as follows. One node's degree increases by one with the probability 1/*N*, because this node is randomly selected. The other node's degree increases by one with the probability *k*_{
i
}/∑_{
j
}*k*_{
j
}, because this node is selected by a random walk from the original randomly-selected node. It is reported that the probability that a walker arrives at this node equals *k*_{
i
}/∑_{
j
}*k*_{
j
}irrespective of the number of steps in the random walk [41]. Note that this probability is equal to that of the probability in preferential attachment [38] which reproduces the heterogeneous connectivity. Thus, the time evolution of *k*_{
i
}is

\frac{\text{d}}{\text{d}t}{k}_{i}=(1-p)\frac{1}{N}+p\left[\frac{1}{N}+\frac{{k}_{i}}{{\displaystyle {\sum}_{j}{k}_{j}}}\right],

(1)

where *N* = (1 - *p*)*t* because the number of nodes increases by one with the probability 1 - *p*, and ∑_{
j
}*k*_{
j
}= 2*t* because one edge is added at every time. Note that this equation is independent of the bypassed path length (the parameter *q*). The solution of the above equation with the initial condition *k*_{
i
}(*t* = *s*) = 1 is

{k}_{i}=[A(p)+1]{\left(\frac{t}{s}\right)}^{p/2}-A(p),

(2)

where *A*(*p*) = 2/[*p*(1 - *p*)].

From the above equation, because *s/t* = *P* (≥ *k*), the cumulative distribution *P* (≥ *k*) is

*P*(≥ *k*) = [*A*(*p*) + 1]^{2/p}[*k* + *A*(*p*)]^{-2/p}.

Since P(k)=-\frac{\text{d}}{\text{d}k}P(\ge k), finally, we get the degree distribution

*P*(*k*) = (*γ* - 1) [*A*(*p*) + 1]^{γ - 1}[*k* + *A*(*p*)]^{-γ},

where the degree exponent *γ* is

\gamma =\frac{2}{p}+1.

(5)

As shown in Equation (4), the degree distribution follows a power law with a cutoff within a small degree.

#### Degree-dependent clustering coefficient

Next, we show an analytical solution for the degree-dependent clustering coefficient of the model via mean-field analysis based on [39].

The clustering coefficient [10, 12] of node *i* is defined as

{C}_{i}=\frac{2{M}_{i}}{{k}_{i}\left({k}_{i}-1\right)},

(6)

where *M*_{
i
}is the number of edges among neighbors of node *i*. Here we consider the time evolution of *M*_{
i
}. The number of edges *M*_{
i
}increases with the probability *p* × *q*, because *M*_{
i
}increases when Event II occurs and a path of length 2 is bypassed (a triangle is generated). That is, we do not need to consider a bypassed path of length greater than 2. Then, *M*_{
i
}of each node, which belongs to the triangle, approximately increases by one. *M*_{
i
}of one node increases by one with the probability 1/*N*, because this node is selected at random. *M*_{
i
}s of the other two nodes increase by one with the probability *k*_{
i
}/∑_{
j
}*k*_{
j
}, because these nodes are selected by a random walk. Therefore, the time evolution of *M*_{
i
}is

\frac{\text{d}}{\text{d}t}{M}_{i}\simeq pq\left[\frac{1}{N}+2\frac{{k}_{i}}{{\displaystyle {\sum}_{j}{k}_{j}}}\right],

(7)

where *N* = (1 - *p*)*t*, and ∑_{
j
}*k*_{
j
}= 2*t*. Moreover, *k*_{
i
}= [*A*(*p*) + 1](*t/s*)^{p/2}- *A*(*p*) as shown in Equation (2).

The solution of the above equation with the initial condition *M*_{
i
}(*t* = *s*) = 0 is

{M}_{i}=2q[A(p)+1]{\left(\frac{t}{s}\right)}^{p/2}+\frac{q(p-2)}{1-p}\mathrm{ln}\frac{t}{s}-2q[A(p)+1],

(8)

where *A*(*p*) = 2/[*p*(1 - *p*)]. From Equation (2), since *k*_{
i
}= [*A*(*p*) + 1](*t/s*)^{p/2}- *A*(*p*), this equation is rewritten as

{M}_{i}=2q\left[\left({k}_{i}-1\right)-\frac{1}{2}\left(2-p\right)A\left(p\right)\mathrm{ln}\frac{{k}_{i}+A\left(p\right)}{1+A\left(p\right)}\right].

(9)

Substituting this equation into Equation (6), we finally get the degree-dependent clustering coefficient

C(k)=q\left[\frac{4}{k}-\frac{2\left(2-p\right)A\left(p\right)}{k\left(k-1\right)}\mathrm{ln}\frac{k+A\left(p\right)}{1+A\left(p\right)}\right].

(10)

#### Average clustering coefficient

Finally, we show a mathematical solution for the average clustering coefficient of the model. Since the average clustering coefficient is expressed as the summation of the product of the degree distribution and the degree-dependent clustering coefficient, it can be described as

C={\displaystyle {\int}_{2}^{{K}_{m}}P\left(k\right)}\times C\left(k\right)\text{d}k,

(11)

where *K*_{
m
}is the maximum degree. The maximum degree is the case that the cumulative probability equals 1/*N* ; thus *P*(≥ *K*_{
m
}) = 1/*N*, and from Equation (3), *K*_{
m
}can be expressed as

*K*_{
m
}= *N*^{p/2}[*A*(*p* + 1)] - *A*(*p*).

Equation (11) is solved via numerical integral because it is analytically unsolvable.

### Estimation of model parameters

This model has two parameters *p* and *q*. In order to reproduce structural properties in metabolic networks, we need to estimate these parameters in real-world networks. In this section, we show how to estimate the parameters.

#### The case of the parameter *p*

Here, we consider the average degree ⟨*k*⟩ of this model.

The average degree is defined as \u3008k\u3009=\frac{1}{N}{\displaystyle {\sum}_{i}{k}_{i}}. As shown in the previous section, *N =* (1 - *p*)*t*, and ∑_{
i
}*k*_{
i
}= 2*t*. That is, the average degree of this model is

\u3008k\u3009=\frac{2t}{\left(1-p\right)t}=\frac{2}{1-p}.

(13)

From this equation, therefore, the parameter *p* is estimated by

p=1-\frac{2}{\u3008k\u3009},

(14)

where ⟨*k*⟩ is obtained from real metabolic networks.

#### The case of the parameter *q*

Here, we consider the number of triangles *T* of this model.

In this model, the number of triangles approximately increases by one with the probability *p* × *q* because a triangle is generated with the probability *q* when Event II occurs. That is,

*T* ≃ *pqt*.

Since *N* = (1 - *p*)*t*, this equation is rewritten as

T\simeq pq\frac{N}{1-p}.

(16)

From this equation, therefore, the parameter *q* is estimated by

q\simeq \frac{T}{N}\frac{1-p}{p},

(17)

where *T* and *N* are obtained from real metabolic networks.

### Data set

We used the metabolic networks of 113 organisms, which were previously investigated in Reference [18]. These metabolic networks are represented by undirected graphs in which nodes and edges correspond to metabolites and substrate-product relationships, respectively. For example, we consider a reaction S1+S2 → P1+P2. In this case, metabolites S1 and S2 connect to products P1 and P2, respectively. That is, the edge list is as follows: (S1, P1), (S1, P2), (S2, P1), (S2, P2). Note that if there are stoichiometric coefficients in the metabolic data used, then they are neglected. In order to accentuate constitutive pathways, these networks exclude 13 ubiquitous metabolites that serve for energy exchange, exchange of a proton or a phosphate moiety, and so on. To be exact, the following metabolites are excluded: water, ATP, ADP, NAD, NADH, NADPH, carbon dioxide, ammonia, sulfate, thioredoxin, (ortho) phosphate (P), pyrophosphate (PP), and H^{+}. We only focused on the largest components of the metabolic networks in order to more accurately evaluate the structural properties.

### Maximum likelihood method considering a cutoff

In order to obtain the degree exponent from real metabolic networks, we used the maximum likelihood method [27]. However, this original method does not consider a cutoff, which we denote by the constant *A*(*p*), in the degree distribution. Thus, it is difficult to compare of the degree exponent between the model and the real data. Consequently, we consider an extended maximum likelihood method:

\gamma =1+N{\left[{\displaystyle \sum _{i=1}^{N}\mathrm{ln}\frac{{k}_{i}+A\left(p\right)}{{k}_{min}+A\left(p\right)}}\right]}^{-1},

(18)

where *k*_{
min
}is the minimum degree in a network.

### Null model

We used a null model to validate our model. The null model is an uncorrelated random scale-free network [28, 29], and is a popular model. Assuming a power-law degree distribution, in the null model, we can obtain a null hypothesis for the degree-dependent clustering coefficient *C*(*k*) and the average clustering coefficient *C* using

C\left(k\right)=C=\frac{{(\u3008{k}^{2}\u3009-\u3008k\u3009)}^{2}}{N{\u3008k\u3009}^{3}},

(19)

where ⟨⋯⟩ denotes the average over all nodes. The values, ⟨*k*⟩, ⟨*k*^{2}⟩, and *N*, are obtained from real metabolic networks.

### Indices for cyclic property

In order to characterize cyclic properties of networks, we define two indices inspired by the cyclic coefficient [30].

One is the cycle index of node *i*, defied as

{r}_{i}^{c}=\frac{2}{{k}_{i}\left({k}_{i}-1\right)}{\displaystyle \sum _{\u3008jh\u3009}{R}_{jh}^{i}},

(20)

where

{R}_{jh}^{i}=\{\begin{array}{ll}1\hfill & (\text{Ifthereisatleastonecycle}\hfill \\ \text{thatpassesthroughnode}i\text{anditstwoneighbors}j\text{and}h)\hfill \\ 0\hfill & (\text{Otherwise}),\hfill \end{array}

(21)

and ⟨*jh*⟩ denotes all pairs of neighbors of node *i*. In addition, *k*_{
i
}is the degree of node *i*. We can understand that this index is an extended clustering coefficient. This index considers cycles whose length is at least 3; however, the original clustering coefficient only focuses on cycles of length 3.

The second index is the cycle length index of node *i*, defined as

{r}_{i}^{l}=\frac{{\displaystyle {\sum}_{\u3008jh\u3009}{({L}_{jh}^{i}-2)}^{-1}}}{{\displaystyle {\sum}_{\u3008jh\u3009}{R}_{jh}^{i}}},

(22)

where {L}_{jh}^{i} is the length of the smallest cycle that passes through node *i* and its two neighbors *j* and *h*.

In order to characterize global cyclic properties, in this section, we focus on the average indices \u3008{r}^{c}\u3009=\frac{1}{N}{\displaystyle {\sum}_{i=1}^{N}{r}_{i}^{c}} and \u3008{r}^{l}\u3009=\frac{1}{N}{\displaystyle {\sum}_{i=1}^{N}{r}_{i}^{l}}, where *N* is the total number of nodes. Small values of ⟨*r*^{c}⟩ indicate a low frequency of cycles in networks. Moreover, small ⟨*r*^{l}⟩ means that the cycle length is globally long in networks.

### Ignoring cycles generated by the network representation

Using these indices as cyclic properties, we investigated the resulting cyclic properties in the metabolic networks of 113 organisms in order to test the hypotheses for cycles, which are generated due to the emergence of the short-cut path. However, we cannot directly use the metabolic networks, which are analyzed in Reference [18] because the metabolic networks include cycles, which are drawn by the network representation. For example, we consider a reaction S1+S2→P1+P2. In this case, a cycle of length 4, is generated as shown in Figure 11.

In this manner, cycles due to the network representation would be drawn when the types of all metabolites are different in a reaction and, as a result, the right-hand side and the left-hand side concurrently consist of multiple metabolites. Therefore, we ignored such cycles when calculated the cycle indices.

### Statistical analysis

In order to assess the significance of the observed correlations, we used Pearson's correlation coefficient *r*, Spearman's rank correlation coefficient *r*_{
s
}, and their *P* -value *P*. We determine that there is a significant correlation between a structural property and optimal growth temperature when *P* < 0.05.