Fingerprint-based method [4]

A chemical compound is represented by a fingerprint of the molecular structure, and the Tanimoto coefficient between fingerprints for compounds is calculated to indicate similarity. It then predicts that there is a metabolic pathway between compounds if the similarity exceeds a certain threshold. The necessary calculations are fast, but accurate path prediction is difficult.

Maximum common substructure search method [5, 6]

This approach focuses on the maximum common substructure between compounds to predict a metabolic pathway. The maximum common substructure search is an NP-hard problem, and requires enormous computation time in order to evaluate the similarity between compounds of complex structures [7, 8]. Various approximation algorithms have been studied in the search for a computationally tractable approach [9–16].

Rule-based method [17–24]

This requires a database of reaction rules constructed from known metabolic reactions, and attempts to predict a metabolic pathway as a sequence of reaction rules. As a feature of reaction rules, some techniques focus on physicochemical properties and structures [25], while other methods focus on enzyme and gene information [26, 27]. Since prediction ability depends on the size and type of metabolites used to build reaction rules, comprehensive prediction is difficult using the exhaustive search algorithm such as breadth-first search [23, 24], and the approach has only been used to predict specific pathways and enzymes. In addition, the complexity of the features used to construct the reaction rules is a factor that has made comprehensive prediction difficult.

This study aims at a comprehensive and de novo approach to predict metabolic pathways between two arbitrary known or unknown compounds, and belongs to the rule-based methods. Using a simple feature that focuses only on the structural formula of compounds, our method enables comprehensive prediction that has been difficult for the conventional methods. Enzyme reactions on the metabolic pathways are used as the reaction rules, and are extracted from the KEGG database [28]. The feature on which we focus is the information before and after the structural change of the compound caused by the enzyme reaction. By using the information about this structural change, our method predicts the enzyme reactions that give the shortest path between two compounds for a given query. Further, a constraint for applying an enzyme reaction rule to a compound is set as the substrate inclusion condition, that is, the compound must include the substrate of the enzyme reaction as part of its own structure. This constraint and shortest-path strategy lead to de novo prediction of unknown biosynthetic pathways that a knowledge-based approach [17, 18] cannot predict. In this study, the metabolic pathway prediction problem is reduced to the shortest path problem, and the search method is based on A* algorithm to traverse nodes in the order of priority and employs the LP solution as an admissible heuristics for estimating the distance to the goal.

KEGG reaction data

The data for compounds and metabolic enzyme reaction information used in this method all come from KEGG. First, we extracted the information pathways from KEGG pathway [29]. By using the KEGG API, we concretely collected all enzyme reactions registered in the global map on the KEGG pathway. In the KEGG Reaction, a pair of compounds that are registered as "main" before and after the reaction indicate that it is a metabolic reaction present in KEGG enzyme or the KEGG pathway global map. In this study, the 2D-SDF structure was extracted only for those pairs registered as "main" to ensure the focus on compound metabolic reactions. Further, most KEGG reactions are registered as reversible reaction, and therefore, the forward and reverse directions were treated as a separate reaction. As a result, the 14570 enzyme reactions and the 6073 related compounds were obtained.

Representation of chemical compounds and enzyme reactions

A key idea in our method is that a chemical compound is converted to a feature vector that represents substructure statistics extracted from the structural formula of the compound. This feature-vector representation evaluates whether a feature, such as a specific substructure, exists in a chemical compound or how many times that feature appears. This converts information about compounds into numerical vectors, called feature vectors, whose i th value corresponds to the existence or frequency of the i th feature considered. This feature-vector representation enables us to reduce the pathway search problem to a computationally feasible problem in the vector space, as will be discussed later in detail.

where ${\mathcal{P}}_l^u$ is a set of paths whose length (depth), or number of bonds, is between l and u (u ≥ l) and which appear at least once in the chemical structures in the dataset. f _c (p) is the number of appearances of path p in the structure of a chemical compound c.

We call the path length range specified by l and u the "representation-depth", and denote it by "depth l-u", which is crucial for the expressiveness of vector representation.

According to the feature-vector representation of chemical compounds, every enzyme reaction rule in the 14570 KEGG enzyme reactions is represented as an operator vector. An operator vector expresses the change in chemical structure before and after the reaction, which is computed as the subtraction of the substrate compound vector from the product compound vector: Let i denote the substrate compound of an enzyme reaction a and j denote the product compound of a. Let $D_l^u\left(i\right)$ and $D_l^u\left(j\right)$ denote the vector representation of i and j, respectively. Then, the operator vector O _a for the reaction a is defined as (see also Figure 1):

Further, every reaction rule R _a for the reaction a is defined as a pair R _a = (U _a , O _a ) of the substrate vector $U_a\left(=D_l^u\left(i\right)\right)$ and the operator vector O _a . As a result, an application of the enzyme reaction to a compound can be achieved simply by "addition" of the operator vector to the compound vector (see also Figure 1), and a reaction pathway from a start compound S to a goal compound G is represented by a sequence of applications (additions) of operator vectors:

$S\underset{U_1\le S}{\overset{R_1=\left(U_1,O_1\right)}{\to }}{X}_1\left(=S+O_1\right)\underset{U_2\le X_1}{\overset{R_2=\left(U_2,O_2\right)}{\to }}{X}_2\left(=X_1+O_2\right)\to \cdots \to X_{n-1}\underset{U_n\le X_{n-1}}{\overset{R_n=\left(U_n,O_n\right)}{\to }}G.$

In this method, different reactions may sometimes be represented by the same vector because of insufficient short-length path counts in the compound vector.

Two constraint conditions for applying enzyme reaction rules

As a constraint for applying a reaction rule to a compound, the substrate inclusion condition is set as inclusion of the substrate vector. When attempting to apply the operator vector of a reaction rule R = (U, O) to a compound vector C, the compound vector C must include the substrate vector U as part of its own structure (see Figure 2). The precise determination procedure requires the entries c _k in the compound vector C = (c₁, ..., c _n ) to exceed the corresponding value u _k in the substrate vector U = (u₁, ..., u _n ), as represented by the following formula:

$\begin{array}{cc}\hfill c_k\ge u_k,\hfill & \hfill \left(1\le k\le n,C=\left(c_1,\dots ,c_n\right),U=\left(u_1,\dots ,u_n\right)\right).\hfill \end{array}$

(4)

Note that this computationally easy procedure for substrate inclusion is a great advantage of our method using vector representation, because the graph inclusion problem for determining whether a compound structure contains a substrate structure is computationally hard (NP-hard).

The second constraint is the "non-negative" compound-vector condition. Since the operator vector O denotes a change in structure before and after the reaction, it sometimes contains negative values and the application (addition) of this to a compound vector C may produce a vector containing negative values. A compound vector is a vector that we define as representing the frequency of occurrence of partial structures, so negative values are not appropriate. Therefore, we set the following non-negative condition to filter out compound vectors containing negative values, preventing inappropriate products that contain negative values as a result of applying an operator vector.

$\begin{array}{cc}\hfill c_k+o_k\ge 0,\hfill & \hfill \left(1\le k\le n,C=\left(c_1,\dots ,c_n\right),O=\left(o_1,\dots ,o_n\right)\right).\hfill \end{array}$

(5)

Search algorithm between two compounds

The purpose of this study is, given a start compound S and a goal compound G, to find the pathway leading to G by applying reaction rules. Here, the distance of a pathway between S and G is defined as "the number of rule applications (i.e., path length)". By introducing this distance definition, the metabolic pathway prediction problem between compounds can be replaced by a mathematical shortest path search problem. That is, finding the shortest path for reaching the integer vector of the goal compound by successively adding integer vectors of reaction rules to integer vectors of intermediate compounds can be considered as the problem of finding the shortest path in the integer-vector space. In this search space, a node is an intermediate compound vector and an edge is an applied reaction operator vector (see Figure 3).

A* algorithm and heuristics

where h(t) is the true distance to the goal. That is, a heuristic function h'(t) that always underestimates the distance to the goal is required. Such a heuristic function h'(t) is referred to as an "admissible heuristics". If a given heuristic is admissible, the A* algorithm will reliably find a shortest path. The A* algorithm was implemented using the data structure "sorted priority queue" for maintaining the nodes to be traversed with weights of the evaluation function value f(t), while the breadth-first search uses the simple "queue" for the nodes to be traversed with no weight.

Datasets and target pathways

DDT degradation pathway

In this study, we used the well-known DDT degradation pathway data set [34] as pathway data to verify the validity of our method. DDT is a chemical substance that can be synthesized for minimal cost, and began to be used as an insecticide during the 1940s because of its insecticidal action against many insects. However, the human carcinogenicity of DDT and its long-term persistency in the environment has since been pointed out [35]. It is important to evaluate the negative impact on the environment, and human health studies analyzing the metabolism of DDT has continued in recent years [36].

Taking into account the number of involved pathways and compounds, as well as the fact that the pathway is a closed circuit, we consider the DDT degradation pathway ideal for verifying our approach. The pathway consists of 20 compounds and 46 enzyme reactions (Figure 4).

Table 2 shows the number of reactions represented as operator vectors associated with the 46 enzyme reactions from the DDT pathway.

Representation-depth	0-1	0-2	0-3
Number of operator vectors	38	44	46

This table shows the number of reaction rules focusing only on the DDT degradation pathway in the KEGG database.

In our experiments, 20 × 19 = 380 pathway routes were selected for the search problem. The first validation experiment only used the 46 enzyme reaction rules contained in the DDT degradation pathway. In the second "more general" experiment, all KEGG reaction rules were used to search the DDT pathway.

Reconstruction of DDT pathway by shortest path finding

Computational times for heuristics

Comparing the efficiency of the heuristic functions in this table showed in particular that a significant reduction in computational time was achieved by the LP heuristic. On the other hand, in the depth 0-3, reduction in computation time was not seen for most heuristics. This implies that, as the representation depth increases, the substrate inclusion condition works more effectively, and the number of branches in the search space becomes smaller.

Prediction of DDT pathway using all KEGG reaction rules

The agreement rate between the true distance and the true pathway route using the LP heuristic were 100% (380/380). Thus, despite using the generic operators (all KEGG reaction rules), the results showed that the method had high reproducibility.

Prediction of Lutein biosynthesis pathway using all KEGG reaction rules

Another pathway prediction using all KEGG reaction rules was executed for Lutein biosynthesis pathway. Lutein biosynthesis pathway is a secondary metabolic pathway from the start compound "Lycopene" to the goal compound "Lutein". Lycopene is a red carotenoid and Lutein is a plant carotenoid, and there are two routes from Lycopene to Lutein in KEGG pathway database: the one is via Zeinoxanthin and the other is via α-Cryptoxanthin. The Lutein biosynthesis pathway has other difficulty compared with DDT pathway prediction: the structures of chemical compounds in the pathway are significantly larger than the ones in DDT pathway, and the KEGG pathway predictive tool PathPred [23, 37] could not predict this pathway.

Our method with the LP heuristics succeeded to precisely predict all pathways between every pair of compounds on the Lutein biosynthesis pathway. The average computational time for the LP heuristic to predict the shortest paths for all pairs was 10.9 seconds. On the other hand, PathPred failed to predict the pathway between Lycopene and Lutein, where the default parameters of PathPred were used: "Simcomp Threshold" was set at 0.4, "Prediction cycle" was set at 1, and Reference pathway was set at "Biosynthesis of Secondary Metabolites (Plants)".

Finding novel biochemical pathways for secondary metabolites of plant origin

To demonstrate the effectiveness of our method for finding novel pathways, we applied our method to predict a biochemical pathway for the start node "Delphinidin" and the goal node "Gentiodelphin". Both compounds are present in the KEGG database. Gentiodelphin is a plant-derived secondary metabolite associated with blue dye, and is known to be synthesized from Delphinidin [23]. The KEGG pathway predictive tool PathPred was also used for performance comparison.

Our method with the LP heuristics predicted the two shortest path solutions shown in Figure 5. Each arrow indicates the enzyme reaction as routing information, accompanied by the KEGG reaction number. Both predicted pathways consist of four enzyme reactions. The first path (blue) is a metabolic pathway present in the KEGG pathway database. The second path (orange) is new and not registered in the KEGG database, and there is the possibility of a new route where the enzyme reaction "R6798" is applied at the end. The computational times for the LP heuristic to predict the shortest paths were 35.0 seconds for the first solution, and 58.4 seconds for the second solution. PathPred only predicted the first pathway with a computational time of 1462 seconds, where the default parameters of PathPred were used.

Overall, our A*-based algorithm with the LP heuristic is more comprehensive and computationally efficient prediction method for biochemical pathway finding.

Search space for pathway predictions

where P is the size of the search space if all solutions are explored, N is the number of reaction rules, and d is the distance from the start node to the goal node. In the heuristic search, the search space can be reduced by visiting the nodes on the true path on a priority basis. Table 5 shows that the LP heuristic can significantly reduce the search space compared to searching all possible solutions.

where B is the ratio that bounds the branching, that is, the ratio at which operator applications are eliminated. When B increases, the base of the exponential function becomes smaller and hence the exponential increase can be reduced. That is, B plays a role in minimizing the exponential expansion of the search space. The significant reduction in computational time achieved by increasing the representation-depth for the vector representation is considered to be due to this reason. In other words, designing a high-speed searching method requires both an accurate heuristic function that estimates the distance to the goal and an effective bound on the branching to reduce the search space.

Reproducibility of KEGG Pathway

Our experimental results for comprehensive predictions using all 8108 KEGG reaction rules show that our proposed method is able to reproduce enzyme reaction pathways in the KEGG pathway database with high accuracy. This is presumably due to the LP heuristic and bound on branching due to the substrate inclusion constraint on the vector representation.

De novo prediction of known and unknown biosynthetic pathways

Our proposed method in this paper is a de novo prediction method in the sense that it does not require any initial network such as KEGG metabolic network as input and it is not a method just to traverse the pathway network. Our method takes as input the set of enzyme reaction rules collected from the KEGG pathway database. However, this does not necessarily imply that the pathway prediction using the list of all reaction rules is equal to the path search on KEGG pathway network. For each compound occurring at a node in KEGG pathway network, the KEGG network only contains the enzyme reactions whose substrate is exactly equal to the compound as an edge connected to the node. On the other hand, our method applies all reaction rules to a given compound if the compound is not only equal to the substrate of the reaction rule but also contains the substrate as a sub-structure (the substrate inclusion condition). Therefore, the search space of our method is exponentially larger than the KEGG pathway network. Further, our method is able to predict unknown biosynthetic pathways between two arbitrary known or unknown compounds.