We choose active main protease (Mpro) and active papain-like protease (PLpro) inhibitors whose pIC50 value are higher than the “activity threshold” as the “seed” set. Eventually, a total of 60 of them are selected, which are denoted as s 1 , s 2 , … , s 60 to be an example. (The active inhibitors are obtained from the articles of Amin et al., 2021 and Ghosh et al., 2021). The inhibitors' molecular structures are represented by SMILES and are shown in Table 1 (only a part of the inhibitors are displayed; all inhibitors' structures with SMILES notations are shown in the Supplementary Material).

There are many measures which can calculate the distance between pairwise inhibitors based their SMILES representation (Weininger, 1988). Most of them use descriptors to extract features and calculate the distance by using classic distance, such as, Manhattan distance, Euclidean distance, Chebyshev distance, and cosine distance. But the design of the descriptor in the feature extraction is not so easy, and it loses some of the information, which we are unsure is useful. On the one hand, some statistical characteristics, such as “SMILES atoms” S _k , the combinations of two “SMILES atoms” SS _k , and the combinations of three “SMILES atoms” SSS _k , take into account the information on the lower-order neighbors of each atom in the molecule, such as the first-order neighbor, second-order neighbor, and third-order neighbor but lack information on the higher-order neighbors of the atom. Moreover, some can also define more optimal descriptors (Toropov et al., 2011), such as BOND, NOSP, and HALO. But those manually designed features only describe part of those inhibitor information, and some unknown important information may still be missed due to the complexity of the feature engineering.

On the other hand, the SARS-CoV PLpro and SARS-CoV Mpro inhibitors can be represented as graphs through their three-dimensional molecular structure. These graphs contain all the topological information of the inhibitor molecules. For these graphs, some metrics, such as the count version of ECFP4 fingerprints, can be used to measure the distance between pairwise inhibitors with different dimensions by extracting features from the graphs. However, these features are not handy for generating new molecules without a structure yet from a set of similar inhibitors. Since we would like to synthesize a novel drug by recombining the structure of a set of highly related molecules, a novel matrix norm was proposed to measure the distance between the pairwise inhibitors with different dimensions without extracting their features. Hence, in this study, a graph representation is conducted to represent a given inhibitor by using its weighted adjacency matrix. The weight is determined by the type of the chemical bond, where the single bond is 1 and the double bond is 2. It is noted that different inhibitors may result in adjacency matrices with different sizes. A graph representation is shown in Figure 3.

Once the pairwise distances between any two inhibitors are obtained by D _p (A,B) ⁱ , i ∈ {1, 2}, a clustering procedure can be conducted to group similar inhibitors, where the shorter the distance between two inhibitors, the higher the possibility that they are grouped into the same cluster. However, not every clustering method works in this case of non-equidimensional data clustering. If the dimensions of two inhibitors are the same, then the cluster center can be naturally obtained in an averagely weighted manner. But if the dimensions of two or more weighted adjacency matrices are different, the center of a group of inhibitors is unavailable by using the above averagely weighted manner. This means, we cannot use the clustering methods that are based on the centroid linkage or rely on the cluster center, such as k-means. In this study, the hierarchical clustering method that D'Andrade (1978) based on the average linkage of two clusters d a v g C i , C j = 1 | C i | ∗ | C j | ∑ x ∈ C i ∑ z ∈ C j d i s t x , z . was employed to test our proposed method.

After clustering, a list of clusters C 1 , C 2 , … , C m can be obtained, and the inhibitors in the same cluster C i = s i 1 , s i 2 , … , s i n can be hybridized to generated new predictions by iteratively separating and reorganizing.

We use a “bottom–up” aggregation strategy to design an iterative algorithm with the heuristic measure function f, which is constructed by D _p (A,B) ⁱ . First, for cluster C = s 1 , s 2 , … , s n , C ∈ C 1 , C 2 , … , C m , each inhibitor s _i is regarded as an initial sample, and then the two closest samples s j ∗ and s k ∗ are found and merged in each step of the algorithm operation.

Then through merging to the hybrid, the closest two inhibitors s j ∗ and s k ∗ in the same cluster by separating and reorganizing, we will get much newer inhibitors s′. Separating means randomly disconnecting a cut edge in the two inhibitor molecules' graphs, which will divide the two inhibitor molecules into four parts. Reorganizing means randomly choosing two of the four parts which do not belong to the same inhibitor molecule, to combine the two parts and making it a new molecule by adding an edge. The schematic diagram is shown in Figure 4.

Finally, an intermediate product s∗ will be chosen by the heuristic measure function f ( s ′ ) = g ( s j ∗ ) ∗ D p ( s j ∗ , s ′ ) i + g ( s k ∗ ) ∗ D p ( s k ∗ , s ′ ) i = w j ∗ ∗ D p ( s j ∗ , s ′ ) i + w k ∗ ∗ D p ( s k ∗ , s ′ ) i from new inhibitors s′ and taken in the place of these two inhibitors.

But the intermediate product s∗ is not the original two inhibitors after all, we therefore set a weight for each inhibitor W = w i | w i = 1 , i = 1,2 , … , n , and as the number of hybridizations increases, the weight of the corresponding inhibitors will be larger. In this way, when calculating the distance, the inhibitors with more hybridization will have a greater distance than before.

With the iteration of the algorithm, the cluster set will remove the original two inhibitors and add an intermediate product until there is only one inhibitor left in the cluster set. This inhibitor is approximately the cluster center of the cluster set. The pseudo code of the proposed algorithm is described as follows. The code is available and can be downloaded from the Internet at https://www.github.com/HenryHan1997/drug_discover.

We get the weighted adjacency matrix from the active main protease (Mpro) and papain-like protease (PLpro) inhibitors' structure. Then, we use D ₂(A,B)² as the distance between the pairwise inhibitors and use the average linkage d _avg (C _i , C _j ) as the distance between the two clusters to cluster the “seed” set by AGNES hierarchical clustering (Kaufman and Rousseeuw, 2009).

Then, we get a tree-like hierarchical structure of the inhibitors according to the average linkage. The threshold is chosen as 0.0012, since it is the elbow position according to the Figure 5. This indicates the distance between the clusters is as large as possible, while the distance within the clusters is as small as possible. After removing clusters less than two inhibitors, 10 clusters are obtained, which are C 1 , C 2 , … , C 10 . They are marked with different colors in Figure 6.

From the results, it can be seen that basically the inhibitors of the same type are still in the same cluster after clustering, i.e., papain-like protease inhibitors 11 − PL, 56 − PL, 64 − PL are in the same cluster, and the main protease inhibitors 5 − M, 13 − M, 14 − M, 15 − M, 16 − M, 19 − M are in the same cluster. This shows that our proposed distance D _p (A,B) ⁱ based on the Kronecker product (p, m)-norm ⋅ p m ⊗ can indeed measure the similarity between pairwise inhibitors of different dimensions.

We chose one cluster, which contains papain-like protease inhibitors 11 − PL, 56 − PL, 64 − PL as an example and used Algorithm 1 with p = 2 and D ₂(A,B)² to discover new inhibitors and count the number of occurrences. Finally, we selected the three most frequent occurrences for analysis, which are shown in Table 2. The new inhibitors are considered to be valid because their SMILES representation can be successfully parsed by the RDKit.

Algorithm 1

Cluster center generation.

To show that our discovered new inhibitors are approximately the cluster centers, we visualized them in a two-dimensional plane. We used principal component analysis (PCA) and sparse PCA to reduce the dimensionality of the distance matrix, which is calculated from these three new inhibitors, intermediate products, and the original seeds by D ₂(A,B)². The results are shown in Figure 7.

From Figure 7A, we can clearly see that the first and second new inhibitors are probably in the center of the cluster, and the third new inhibitor does not perform well; from Figure 7B, it is evident that the three new inhibitors are probably in the center of the cluster. On the whole, we calculated the sum of the distance from the new inhibitor to the seeds, and showed that the first new inhibitor performs best, which is shown in Figure 8.

Finally, we calculated the quantitative estimate of drug-likeness (QED) of the new inhibitors and original inhibitors, which is synthesized by using eight descriptors. The descriptors contain MW, logP, HBA, HBD, PSA, ROTB, AROM, and ALERTS (Brown et al., 2019). The QED of the first new inhibitor reached 0.731, which is the highest and is higher than the original three inhibitors (11 − PL, 56 − PL, 64 − PL). The results are shown in Table 3.

At the same time, we also record the synthetic route of the first new inhibitor for analysis, which is shown in Figure 9. The first new inhibitor is obtained by recombining papain-like protease inhibitors 56 − PL and 64 − PL to form an intermediate product c1(c(ccc(c1)N)C)COC, and then separating and combining the intermediate product c1(c(ccc(c1)N)C)COC and papain-like protease inhibitor 11 − PL.

This kind of procedure is not possible by using the SMILE–based method directly. But using the proposed graph representation, we can easily generate more number of potential new drugs by combining information of currently known related drugs.