Unveiling the molecular mechanism of SARS-CoV-2 main protease inhibition from 137 crystal structures using algebraic topology and deep learning.

Currently, there is neither effective antiviral drugs nor vaccine for coronavirus disease 2019 (COVID-19) caused by acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Due to its high conservativeness and low similarity with human genes, SARS-CoV-2 main protease (Mpro) is one of the most favorable drug targets. However, the current understanding of the molecular mechanism of Mpro inhibition is limited by the lack of reliable binding affinity ranking and prediction of existing structures of Mpro-inhibitor complexes. This work integrates mathematics (i.e., algebraic topology) and deep learning (MathDL) to provide a reliable ranking of the binding affinities of 137 SARS-CoV-2 Mpro inhibitor structures. We reveal that Gly143 residue in Mpro is the most attractive site to form hydrogen bonds, followed by Glu166, Cys145, and His163. We also identify 71 targeted covalent bonding inhibitors. MathDL was validated on the PDBbind v2016 core set benchmark and a carefully curated SARS-CoV-2 inhibitor dataset to ensure the reliability of the present binding affinity prediction. The present binding affinity ranking, interaction analysis, and fragment decomposition offer a foundation for future drug discovery efforts. This journal is © The Royal Society of Chemistry.

Introduction

Starting in late Dec, 2019, the COVID-19 pandemic caused by new severe acute respiratory syndrome coronavirus (SARS-CoV-2) has infected more than 22 million individuals and has caused more than 777 000 fatalities in all of the continents and over 213 countries and territories by August 19th, 2020. Under the current global health emergency, researchers around the world have engaged in the investigation of the different drug targets of SARS-CoV-2, such as the main protease (Mpro, also called 3CLpro), papain-Like protease (PLpro), RNA-dependent RNA polymerase (RdRp), 5′-to-3′ helicase protein (Nsp13) to seek potential cures for this serious pandemic. To date, although there are some vaccines undergoing the Phase III trials,[1] their safety and efficacy are still unclear.[2]

The main protease, one of the best-characterized targets for coronaviruses, attracts lots of research attention because it is very conservative and distinguished from any human gene. A recent study shows that although the overall sequence identity between SARS-CoV and SARS-CoV-2 is just 80%, the Mpro of SARS-CoV-2 shares 96.08% sequence identity to that of SARS-CoV.[3] Therefore, we hypothesize that a potent SARS Mpro inhibitor is also a potent SARA-CoV-2 Mpro inhibitor.

At this moment, more than 300 potential SARS-CoV Mpro inhibitors with its binding affinities are available in ChEMBL database[4] which can be considered as the potential SARS-CoV-2 Mpro inhibitors. Recently, total 146 crystal structures of SARS-CoV-2 Mpro with its ligand complexes are released on the Protein Data Bank (PDB).[5] Among them, 137 crystal structures have no available binding affinities reported for various reasons. However, the central dogma of drug design and discovery concerns the molecular mechanism and binding affinity of drug target interactions. Knowing the binding affinities and their ranking of 137 SARS-CoV-2 Mpro inhibitors is of great significance to the future design of anti SARS-CoV-2 drugs.

In this work, for the first time, we predict the binding affinities of these 137 Mpro–inhibitor complexes by reformulating algebraic topology-based mathematics-deep learning (MathDL) models, which have been the top competitor in D3R Grand Challenges, a worldwide competition series in computer aided drug design in the past three years.[6] We generate reliable poses for 141 Mpro inhibitors with binding affinities but without complex structures. Together with 44 other complexes, we compose a set of 185 Mpro–inhibitor complexes, which is paired with 17 382 protein–ligand complexes in PDBbind 2019 general set. These datasets are utilized to construct 11 MathDL models in single-task and multitask settings.[6] One of these 11 MathDL models has been validated by using the PDBbind v2016 core set benchmark, achieving the top performance over all exiting scoring functions. The other ten MathDL models have cross-validated on a set of 185 Mpro–nhibitor complexes, showing an averaged Pearson's correlation coefficient of 0.73.

Notably, for covalent inhibitors, the scheme of covalent irreversible inhibition of SARS-CoV/SARS-CoV-2 Mpro is presented below:

The inhibitor first binds to the protease noncovalently, then a nucleophilic attacking by Cys145 leads to the formation of a stable covalent bond between the protease and the inhibitor.[7,8] The interaction depends on both the equilibrium-binding constant Ki (designated as k1/k2) and the inactivation rate constant for covalent bond formation k3. In this work, the binding affinity/IC50 assesses the first step to form noncovalent binding.

In a nutshell, the present work provides reliable binding affinity predictions and ranking of 137 SARS-CoV-2 inhibitors that have crystal structures. It also offers data curation and validated models for exploring potential SARS-CoV-2 Mpro inhibitors. Furthermore, this work explores different possible binding regions on the SARS-CoV-2 main protease and decode the most favorable molecular fragments for the inhibitor design.

Results and discussions

Results

This section is devoted to the utilization of our MathDL models developed in Section 3.3 to predict the binding affinities and their ranking of SARS-CoV-2 inhibitors that do not have reported experimental affinities. To reduce the role of 3D pose prediction errors in our model, we use the SARS-CoV-2 inhibitors with X-ray structures available in the PDB for our study. We manually search these ligands on the PDB and arrive at a set consisting of 137 SARS-CoV-2 Mpro inhibitors having X-ray crystal structures but lacking of experimental binding affinities. We name this set SARS-CoV PDB-noBA (see Table 3). In this experiment, we develop a MathDL model optimized from PDBbind v2016 core set (see Section 3.3.1), five MathDL-ALL and five MathDL-MT models obtained from 5-fold study on the SARS-CoV BA set (see Section 3.3.2). The final predicted binding affinity is the consensus of these 11 models. The top ten inhibitors indicated by our models are shown in Table 1.

A summary of our selected data sets

Data name	Data size	Descriptions	References
PDBbind v2019	17 382	Partial PDBbind general set v2019	18
PDBbind v2016 core set	290	PDBbind v2016 core set	18
SARS-CoV PDB	192	Inhibitors of SARS-CoV/SARS-CoV-2 Mpro having X-ray crystal structures	5, 19 and 20
SARS-CoV PDB-BA	44	Inhibitors of SARS-CoV/SARS-CoV-2 Mpro having X-ray crystal structures and experimental binding affinities	5, 18–23
SARS-CoV PDB-noBA	137	Inhibitors of SARS-CoV-2 Mpro having X-ray crystal structures but lacking of experimental binding affinities	5, 18–20, 24
SARS-CoV 2D	141	Inhibitors of SARS-CoV/SARS-CoV-2 Mpro having only 2D structures	4, 19, 20, 25–31
SARS-CoV BA	185	Inhibitors of SARS-CoV/SARS-CoV-2 Mpro having experimental binding affinities	5, 18–20, 26–31

Binding affinities of top 10 complexes in SARS-CoV PBD-noBA dataset predicted by our MathDL. “Pred. BA” indicates the predicted binding free energy in kcal mol−1 and “Pred. IC50” is the corresponding IC50 in μM unit via the following conversion: Pred. IC50 = 10Pred. BA/1.3633 × 106

PDBID	Pred. BA	Pred. IC50	PDBID	Pred. BA	Pred. IC50
7c8t	−8.90	0.30	6z2e	−8.43	0.66
5rgl	−8.50	0.58	6xbi	−8.34	0.76
6xhm	−8.50	0.58	6xmk	−8.33	0.78
7bqy	−8.49	0.59	5rh7	−8.32	0.79
5rfr	−8.45	0.63	6xbh	−8.27	0.86

The most potent SARS-CoV-2 inhibitor found by our MathDL models is the inhibitor Nol in complex 7c8t. Nol was synthesized by Yang and his colleagues,[9] Nol is found remarkable activities against SARS-CoV and HCoV.[9] Specifically, the dissociation constant Ki of Nol was found to be 0.053 μM against SARS-CoV.[9] Our MathDL reveals that Nol still inhibits SARS-CoV-2 main protease with a potent affinity at 0.30 μM.

Another important top potent SARS-CoV-2 inhibitor found by our models is the Michael acceptor inhibitor N3 in complex 7bqy. Designed by Yang and his colleagues,[8] N3 was found to have viral activities against different coronavirus Mpro such as SARS-CoV and MERS-CoV.[8,10] Specifically, the dissociation constant Ki of N3 was found to be 9.0 μM against SARS-CoV.[8] Our MathDL reveals that N3 still inhibits SARS-CoV-2 main protease with an even better affinity at 0.59 μM. This finding is consistent with the literature work[11] showing that N3 is a potent inhibitor of COVID-19 virus Mpro.

The inhibitor Qys in the complex 6xmk is also noticeable. Our predicted IC50 is 0.78 μM. Soon after we made the prediction, on August 12th, 2020, Rathnayake et al.[12] released another Qys-main protease complex with PDB ID 6w2a and also reported the IC50 of Qys to SARS-CoV-2 is 0.45 μM, which is close to our prediction.

It is worth pointing out, except for the inhibitor T9j in the complex 5rg1, the rest of inhibitors reported in Table 1 are covalent inhibitors, which irreversibly form covalent bonds with Cys145 of the main protease (see discussion in Section 2.2.2). However, our models only predict the non-covalent binding affinity which is measured before the enzyme deactivation. The predicted binding affinities of all 137 complexes in SARS-CoV PBD-noBA dataset from various MathDL models are presented in Table S8 in ESI.† In this table, we also supply the synthetic accessibility score (SAS), partition coefficient log P, and solubility log S for each small molecule. Except for SAS obtained via RDKit,[13] log P and log S are evaluated by our TopP-S model.[14]

Discussion

Binding site analysis

Based on the crystal structure information of 137 complexes in SARS-CoV PDB-noBA set, we have identified 13 distinct binding site regions of the SARS-CoV-2 main protease as illustrated in Fig. 1. Those binding pockets are denoted by P, i = 1, 2, …, 13. Fig. 2a reveals that binding pocket P1 is the most common binding region of the SARS-CoV-2 main protease, which attracts around 80.2% of ligands in the SARS-CoV PDB-noBA data set of 137 complexes. This finding is no surprise since the binding pocket P1 shares similar active sites to its predecessor, i.e. SARS-CoV Mpro. Specifically, P1 encompasses His141 and Cys145 catalytic dyad which are imperative to the substrate-binding mechanism.[8] In additions, the substrate-binding residues Tyr161 and His163 (ref. 15) are covered in P1. Binding pockets P2, P3, P5, P7, P8, and P10 are the least favor sites consisting of only one ligand. The rest of the binding pockets involve no more than 7 ligands. To study the correlation of the binding regions to the binding free energy, we present the box plot in Fig. 2b to illustrate the energy values through their quartiles.

Fig. 1

All binding site pockets observed from 137 inhibitors in SARS-CoV PDB-noBA set.

Fig. 2

(a) Distribution of 137 ligands across 13 distinct binding sites; (b) Box plot of predicted binding energies (kcal mol−1) of all inhibitors in each binding site.

The prevailing binding pocket P1 is the best region on the SARS-CoV-2 Mpro for inhibitor design with the median binding energy being −7.22 kcal mol−1. Nol is the best inhibitor candidate for the binding site P1 with predicted affinity found to be −8.90 kcal mol−1. Other binding regions such as P4, and P11 are less common but show their adequate effects on the binding mechanism with their best energy binding affinities calculated at −7.28 kcal mol−1 and −6.80 kcal mol−1, respectively. These potential binding sites can guide drug combination to inhibit coronavirus Mpro effectively.

Interaction analysis

By looking further into the interactions between the top inhibitors and the main protease, we have found that Nol, V2m, N3 are peptidomimetic inhibitors, they form as many as 8, 8, 9 hydrogen bonds respectively to the nearby residues and also all form 1 covalent bond with Cys145 as listed in Table 2 and depicted in Fig. 3. All these hydrogen bonds justify their potency of the first step of noncovalent binding to the main protease complex and confirms the robustness of our MathDL models; the covalent bonds make the binding irreversible. We also notice that these three inhibitors share two common hydrogen bonds to His163, His164 (see Table 2, Fig. 3a, c and d). Therefore, they have some similar predicted binding energies, especially 6xhm and 7bqy at −8.50 kcal mol−1 and −8.49 kcal mol−1, respectively.

Interaction analysis in the binding pockets of top 4 complexes in term of binding affinity predicted by our MathDL models

PDB ID	Ligand ID	Hydrogen bond	Covalent bond
7c8t	Nol	His163, His164, Cys145, Gln189, Gly143, Glu166	Cys145
5rg1	T9J	His163, Glu166
6xhm	V2m	His163, His164, Cys145, Gln189, Phe140	Cys145
7bqy	N3	His163, His164, Cys145, Gln189, Thr190, Glu166, Phe140, Gly143	Cys145

Fig. 3

The interactions between the top 4 inhibitors in the SARS-CoV PBD-noBA dataset and SARS-CoV-2 Mpro: (a) 7c8t; (b) 5rg1; (c) 6xhm; and (d) 7bqy. Inhibitors are shown in the purple color. Hydrogen bonds are marked in dashed green lines, and covalent bonds are depicted in solid blue lines. All interactions are shown with the distance information in Å.

This examination manifests how well our models preserve and capture the physical and chemical properties described in intermolecular bonding interactions. Furthermore, the ligand T9J that binds to Mpro in complex 5rg1 with a quite close binding energy at −8.50 kcal mol−1 forms different hydrogen bonds in comparison to three previously mentioned inhibitors (see Table 2). Since our models only concern the non-covalent binding affinity, the lack of covalent bond in 5rg1's interactions does not downgrade its binding strength. With two relatively large hydrogen bonding distances (O2-His163: 3.05 Å, O3-Glu166: 3.38 Å (see Fig. 3d)), the binding affinity of 5rg1 is still comparable to the top inhibitors indicating the important roles in acquiring the hydrogen bonds to these residues in the main protease's binding process.

In the top 10 inhibitors as listed in Table 1, T9J in the complex 5rg1 is only one non-covalent inhibitor. The rest belongs to the class of targeted covalent inhibitors (TCI) in which they interacts with the protein residues, i.e., cysteine, to form a covalent complex strongly neutralizing target's function. However, the major disadvantage of TCIs is the association with the high toxicity risks.[16] TCIs' strong covalent bond can irreversibly modify the unintended protein targets in the human body. As a result, the top covalent inhibitors in SARS-CoV PBD-noBA dataset may have little chance to become approved market drugs in comparison to their non-covalent counterparts such as T9J in 5rg1.

Due to the popularity of the binding site P1 among 137 interested inhibitors, we mainly analyze the interaction network around the residues in that region. Out of 110 molecules binding to P1, there are 103 inhibitors forming at least one hydrogen bond to the nearby amino acid in the SARS-CoV-2 main protease. We have identified 20 different residues in the binding pocket P1 composing hydrogen bonds to these small molecules. Fig. 4 illustrates the frequency of these 20 residues across 110 inhibitors. Based on Fig. 4, Gly143 residue is the most attractive site to form the hydrogen bond. It appears in 53.6% of 110 intermolecular bonding interactions, followed by Glu166 residue with a frequency of 39.1%; residue Cys145 and His163 also occupy 38.2% and 30.9%, respectively. It is worth noting when these molecules form a hydrogen bond with Cys145, they also constitute another hydrogen bond with Gly143. In all cases, both these residues share the same hydrogen-bond acceptor. Besides the hydrogen bond network, 71 ligands in the SARS-CoV PDB-noBA dataset form a covalent bond to γ-sulfur of Cys145. Except the second one, all the others in the top 10 inhibitors are equipped with that covalent bond (see Table S8 in ESI†).

Fig. 4

Popularity of amino acids in the binding site P1 constituting the hydrogen bonds with ligands.

Furthermore, we are interested in the binding energy distribution associated with the interaction network. Fig. 5 depicts the violin plot of that distribution across four categories, namely no H-bond (no hydrogen bond), H-bond (at least one hydrogen bond), no cov. bond (no covalent bond), and cov. bond (at least one covalent bond). Hydrogen bond interactions that are expected to play an important role in the binding mechanism are well captured in our MathDL models. Specifically, while the average energy of inhibitors having no hydrogen bond is −6.62 kcal mol−1, the average energy of ones with hydrogen bond is as low as −7.23 kcal mol−1.

Fig. 5

Violin plot of the predicted binding energies for 110 inhibitors binding to the binding site P1 classified into 4 categories, namely no H-bond (no hydrogen bond), H-bond (at least one hydrogen bond), no cov. bond (no covalent bond), cov. bond (at least one covalent bond). The mean is in the orange color, the median is in the blue color, and the minimal and the maximal values are both in the black color.

It is noted that our MathDLs only measure the non-covalent binding affinity. The covalent bond appearing at the final covalent complex is not properly accounted for in our framework. Therefore, it is expected that our models sometimes overestimate the covalent-bond inhibitors over the non-covalent-bond candidates. Fig. 5 reveals molecules in the group of covalent bonds generally are predicted with lower binding energy with an average being −7.42 kcal mol−1 in comparison to −6.89 kcal mol−1 averagely measured on ones without covalent bonds.

Fragment analysis

To design the lead molecules, it is of importance to have promising fragments from existing inhibitors against the drug targets. Therefore, in the present work, we study all the fragments decomposed from 110 inhibitors attached to the binding site P1. To carry out this task, we utilize BRICS algorithm[17]via RDkit.[13] In BRICS model, there are 16 chemical environments indicated by linkers denoted by L1, L2, …, L16. The BRICS decomposition gives raise to a total of 185 unique fragments, which are all presented in Table S9 in ESI.†Fig. 6 illustrates top 12 common fragments in terms of their frequencies. Noting that the second frequent fragment, L1–C(C)O, often constitutes a hydrogen bond with Gly143 and in many cases forms a covalent bond with Cys145.

Fig. 6

Fragment frequencies based on BRICS decomposition of 110 inhibitors of binding site pocket P1. L is the link atom of a certain type described in ref. 17.

Materials and methods

Datasets

Our deep learning-based scoring function, MathDL, was trained on public databases including PDBbind[18] and ChEMBL.[4] The PDBbind sets contain all complexes with crystal structures deposited in the PDB with the binding affinities not limited to Kd, Ki, and IC50 reported in the literature. In this work, we employ the PDBbind v2019, the latest version of its generation. The v2019 version of the PDBbind consists of 17 679 protein–ligand complexes. However, the data preprocessing of the MathDL[32] only retains 17 382 complexes. Among them, there are 10 485 ligands measured in Kd/Ki and 6537 ligands measured in IC50.

ChEMBL is another manually curated database of bioactive molecules. Currently, ChEMBL contains more than 2 million compounds in the SMILES string format. Excluding 30 main protease inhibitors in PDBbind data, we have found other 277 small molecules on ChEMBL with reported Kd/IC50. Additionally, we have found more than 300 other SARS-CoV main protease inhibitors from literatures.[18-20,25-31] In total, there are more than 600 ligands bound to SARS-CoV/SARS-CoV-2 main protease having the experimental binding affinities; among them, there are 44 crystal structures. For compounds without the crystal structures, MathPose[6] is utilized to generate their 3D conformations. The predicted 3D coordinates of these structures are presented in the SDF format and available in ESI.† Currently, there are roughly 137 ligands forming crystal complexes with the SARS-CoV-2 main protease on PDB without the report of the experimental inhibitor activities. Most of them are deposited by the PanDDA analysis group (https://pandda.bitbucket.io/#).

To serve model validation purposes, we classify the selected data into five different groups as listed in Table 3. Specifically, PDBbind v2019 is the biggest set in this compilation with its PDB IDs and experimental binding affinities listed in Table S1 in ESI.† PDBbind v2016 core set is a subset of PDBbind v2019 and is formed by 290 complexes representing all protein classes in the refined set of PDBbind v2016.[18,33] The PDB IDs of all complexes in the PDBbind v2016 core set are provided in Table S2.† We also collect all Mpro complexes of SARS-CoV/SARS-CoV-2 on the PDB, denoted by SARS-CoV PDB, which results in a total of 192 structures (see Table S3†). Among them, there are 44 ligands with the report of experimental binding affinities denoted by SARS-CoV PDB-BA (see Table S4†). Furthermore, we are interested in the set of SARS-CoV-2 Mpro complexes in the aforementioned SARS-CoV PDB set but their affinities are not presented or undisclosed. We call this set SARS-CoV PDB-noBA with PDB IDs listed in Table S5.† To enrich our training data targeting SARS-CoV/SARS-CoV-2 main protease inhibitors, we gather some inhibitors reported on the literature.[4,25] For those compounds with only 2D information, we limit ourselves to ones having the similarity score based on the path-based fingerprint FP2 no lower than 0.6 to at least one inhibitor in the SARS-CoV PDB set. As a result, we arrive at a set of 141 structures named SARS-CoV 2D (see Table S6†). Combining SARS-CoV PDB-BA and SARS-CoV 2D data sets, we finalize a reliable database focusing on SARS-CoV/SARS-CoV-2 main protease inhibitors. Notice that the binding affinities in this set are all reported in IC50. Table S7 in ESI† presents the PDB IDs as well as the experimental binding energies of these ligands.

Methods

MathDL

The MathDL models developed in this work are reformulated from our early model bearing the same name. MathDL was designed for the prediction of various druggable properties of 3D molecules.[6] In the past three years, MathDL has been proved to be the top competitor in D3R Grand Challenges (https://drugdesigndata.org/about/grand-challenge), a worldwide competition in computer-aided drug design. In the present work, we have, for the first time, develop a multitask MathDL (MathDL-MT) to handle the Mpro inhibitor dataset. We have also extended our earlier MathDL by including all different datasets (MathDL-All). Fig. 7 depicts the framework of the MathDL in which the element-specific algebraic topological representations are integrated with the convolutional neural network (CNN) aiming to predict varied druggable properties such as toxicity, binding affinities, etc.

Fig. 7

A framework of MathDL energy prediction model which integrates advanced mathematical representations with sophisticated CNN architectures.

Algebraic topology-based representations

Algebraic topology studies the topological spaces with the use of abstract algebra, which can dramatically simplify the geometric complexity. Persistent homology (PH) is one of the algebraic topology approaches which has the capacity to track the multiscale topological information over different scales along with filtration by characterizing independent components, rings, and higher dimensional voids in space.[34] In this section, we will briefly review the algebraic topology-based representations. Additionally, since we are dealing with the protein–ligand system, therefore, the biological considerations will take into account as well.

Simplex

The q-simplex denoted as σ is the convex hull of q + 1 affinely independent points in . For example, the 0, 1, 2, and 3-simplex is considered as a vertex, an edge, a triangle, and a tetrahedron, respectively. We call the convex hull of each non-empty subset of q + 1 points the face of σ, and each points are also called the vertices.

Simplicial complex

A set of simplices is a simplicial complex denote K which satisfies that every face of a simplex σ ∈ K is also in K and the non-empty intersection of any two simplices in K is the common face for both.

Chain complex

A formal sum of q-simplices in simplicial complex K with coefficients in an algebraic field (typically ) is a q-chain. A set of all q-chains of the simplicial complex K equipped with an algebraic field is called a chain group and denoted as C(K). The boundary operator is defined by ∂: C(K) → C(K) to relate the chain groups. More specifically, we denote σ = [v0, v1, …, v] for the q-simplex spanned by its vertices, and then the boundary operator can be represented as:Here, is the (q − 1)-simplex with vi being omitted. The sequence of chain groups connected by boundary operators is called the chain complex and expressed as:The q-cycle group Z(K) and the q-boundary group B(K) are defined as Z(K) = ker(∂) = {c ∈ C(K)|∂c = ∅} and B(K) = im(∂) = {∂c|c ∈ C(K)}. The q-th homology group is the quotient group H(K) = Z(K)/B(K). Moreover, the rank of q-th homology group can be computed as rankH(K) = rankZ(K) − rankB(K), which is denoted as the q-th Betti number β. To be notice that the q-th Betti number count the number of q-dimensional holes that can not be continuously deformed to each other.

Persistent homology

A filtration of a simplicial complex K is a nested sequence of subcomplexes of K such that ∅ = K0 ⊆ K1 ⊆ K2⋯⊆ K = K. Then the p-persistent qth homology group of K is defined as:Here the rank of H(K) counts the number of q-dimensional holes in K that are still alive in K, which is called the p-persistent qth Betti number. The persistent homology not only records the topological information at a specific configuration, but also tracks the changes along with the filtration parameters. More specifically, the topological changes will be preserved in the persistent barcodes. In MathDL, we make use of the persistent homology barcodes by dividing them into bins and calculating the birth, death, and persistence incidents in each bin to enrich our algebraic topological representations.

Element specific considerations

The protein–ligand complex is structural and also biological. The persistent homology provides a theoretical approach to encode high-dimensional spatial data of protein–ligand complexes into algebraic topological representations. In this section, we address the biological considerations for biomolecular complexity. There are many kinds of interactions that exist in the protein–ligand complex, such as electrostatics, hydrogen bonds, and hydrophobic effects. Although persistent homology can capture the interactions among the nearest neighbors, the long-range interactions will be hindered. This difficulty can be avoided via the deployment of the element-specific attention.[32] There are 4 commonly atom types in protein, namely C, N, O, S, and there are 11 commonly atom types in ligand, including C, N, O, S, P, F, Cl, Br, I, H, B. We include Boron in the ligand atom type consideration since it appears in more than 200 small compounds in our training data. The general framework of MathDL is depicted in Fig. 7 under exemplified steps. In addition, the details of the deep learning architecture of the current MathDL is offered in Fig. S1.† For the details of feature descriptions as well as the deep learning architecture, interested readers are referred to our previous work.[32]

MathPose

MathPose, a 3D pose predictor that converts SMILES strings into 3D poses with references of target molecules, is the top performer in D3R Grand Challenge 4 (GC4) in predicting the poses of 24 beta-secretase 1 (BACE) binders.[6] For one SMILES string, around 1000 3D conformations can be generated by various docking software tools such as GOLD,[35] Autodock Vina,[36] and GLIDE.[37] Moreover, a selected set of known complexes is re-docked by three aforementioned docking software packages to generate at least 100 decoy complexes per input ligand used in the machine learning training set. The machine learning labels will be the calculated root mean squared deviations (RMSDs) between the decoy and native structures for the training data of the pose selection task. Furthermore, MathDL models will be set up and applied to select the top-ranked pose for the given ligand. Besides the GC4 challenge, our models have outperformed state-of-the-art scoring functions at the docking power challenge on CASF-2007 and CASF-2013 benchmarks.[33] Those established results attest to the credibility of our MathPose on the 3D structure prediction of small molecules.

Validations

PDBbind v2016 core set benchmark

In this validation task, we will testify our model against 290 complexes in the PDBbind v2016 core set. This is a prevalent test set to assert the scoring ability of a binding affinity prediction model and has attracted lots of research groups to devote the effort to improve the Pearson's correlation coefficient (Rp) and Kendall's tau (τ) on this core set performance.[18,42,43] In the current work, we merge the PDBbind v2019, SARS-CoV PDB-BA, and SARS-CoV 2D sets but removing the duplicates and excluding the PDBbind v2016 core set complexes to attain a training set of 17 211 complexes. MathDL with the architecture described in Section 3.2.1 is trained on those complexes. The resulting model is utilized to predict the binding affinity of 290 structures in the PDBbind v2016 core set.

With the purpose of exploring the most optimal model for this benchmark, MathDL is trained for 1000 epochs. Then, we pick the epoch based on the root-mean-squared error (RMSE) of the PDBbind v2016 core set prediction. We have found that MathDL achieves the smallest RMSE in this experiment at 140 epochs. Specifically RMSE, Rp, and τ metrics on the v2016 core set are 1.56 kcal mol−1, 0.858, and 0.671, respectively. Meanwhile, the training accuracy is 0.387 kcal mol−1 in terms of RMSE and its Pearson's correlation coefficient is Rp = 0.994. These performances reveal that our MathDL converges very fast and with only 140 epochs and maintains a good balance between training and testing accuracies. This is a state-of-the-art performance since our MathDL is ranked in the second place in comparison to 33 other scoring functions (see Fig. 8). It is noted that the top model is TopBPcon. published in our previous work[32] with Rp = 0.861. TopBPcon. is the consensus of gradient boosted tree and deep learning-based models. If only the deep learning framework is considered, the performance of TopBP (denoted by TopBP-DL) on the core set of PDBbind v2016 is Rp = 0.848.

Fig. 8

The Pearson correlation coefficient of various scoring functions on PDBbind v2016 core set benchmark. Our MathDL is in red. The performances of other models that are in green are taken from ref. 18, 32, 38–41. TopBPcon., the consensus model in our published work,[32] attains the highest Rp at 0.861. The current MathDL is followed with the second highest Rp at 0.858 and RMSE = 1.56 kcal mol−1. The third place in the list is another TopBP model, TopBP-DL, solely based on the deep learning architectures and its reported Rp is 0.848.[32] It is noted that all of the machine learning based scoring functions in this comparison were trained on the PDBbind v2016 refined set of 3767 complexes except for our MathDL. Explicitly, MathDL is trained on a much larger training set consisting of 17 211 complexes picked out from the PDBbind v2019 set and SARS-CoV BA set.

It is worth mentioning that except for our MathDL, all machine learning-based scoring functions listed in Fig. 8 were trained on the PDBbind v2016 refined set of 3767 complexes. As mentioned above, the current MathDL is compiled on a much larger training set comprised of 17 211 complexes selected from PDBbind v2019 and SARS-CoV BA data. Even the present MathDL has not outperformed its predecessor, i.e., TopBPcon., MathDL is still a preference model since it is trained on a diverse data set covering various protein families and different binding energy ranges. As a result, it is expected to deliver more reliable predictions on the SARS-CoV-2 inhibitor, especially when this main protease family is not included in the training data of previous TopDL models. The resulting MathDL model is labeled as MathDL-Core2016 and is utilized to predict affinities of complexes in SARS-CoV PDB-noBA in Section 2.1.

5 fold cross-validation on SARS-CoV BA set

In this section, we testify the performance of our MathDL against 185 inhibitors in the SARS-CoV BA set aforementioned in Table 3. Among those ligands, there are 44 X-ray crystal structures and the rest are in 2D SMILES strings. We employ MathPose to predict 3D structures of those 2D ligands. To carry out the validation, we randomly split the SARS-CoV BA set into 5 non-overlapped folds. In each fold prediction task, MathDL trains on the partial data of SARS-CoV BA in conjunction with PDBbind v2019 set. This situation results in two different ways of training our MathDL model. The first approach is a traditional MathDL architecture with the training set combining both SARS-CoV BA and PDBbind v2019 complexes. The second model makes use of multi-task learning.[44] In each epoch, the weights of the MathDL architecture are learned through the information from PDBbind v2019 set, then only the fully connected layers are trainable when learning SARS-CoV BA structures. Finally, we come up with 10 different MathDL models in which the traditional MathDL frameworks are labeled as MathDL-All-i and multi-task MatDL is named MathDL-MT-i with i running from 1 to 5. In each model, after 100 epochs, we start monitoring which epoch that helps our model achieve the smallest RMSE on the test set.

Table 4 reveals that MathDL-All models are well trained with the averaged accuracy RMSE = 0.286 kcal mol−1, Pearson's correlation coefficient Rp = 0.994, and Kendall's tau τ = 0.934. Their averaged performances on test data across 5-fold of the SARS-CoV BA set are found to be Rp = 0.729, τ = 0.540, and RMSE = 0.789 kcal mol−1. These results endorse the reliability of these models in the binding affinity prediction of SARS-CoV/SARS-CoV-2 inhibitors. Table 4 also lists the training and testing performances of five multi-task learning models. The averaged training performance of the MathDL-MT model is Rp = 0.995, τ = 0.941 and RMSE = 0.275 kcal mol−1. The accuracy of the multi-task architecture on the test sets is similar to MathDL-All with Rp = 0.727, τ = 0.532, and RMSE = 0.822 kcal mol−1. With these promising results, it is encouraging to carry out MathDL models to predict unknown binding affinities of SARS-CoV/SARS-CoV-2 inhibitors. It is worth noting that if the 5-fold cross-validation is conducted purely on the SARS-CoV BA set, the average Rp and τ are as low as 0.561 and 0.388, respectively. These results strongly support the inclusion of diverse information such as PDBbind v2019 in conjunction with sophisticated deep learning architectures to achieve the accurate binding energy prediction of Mpro inhibitors.

5-fold Performances of MathDL-All and MathDL-MT on SARS-CoV BA set

	MathDL-ALL			MathDL-MT
	R p	τ	RMSE	R p	τ	RMSE
Fold 1 (train)	0.992	0.923	0.327	0.996	0.949	0.253
Fold 1 (test)	0.792	0.534	0.682	0.818	0.534	0.680
Fold 2 (train)	0.995	0.943	0.266	0.996	9.948	0.236
Fold 2 (test)	0.625	0.498	0.866	0.689	0.538	0.826
Fold 3 (train)	0.991	0.917	0.367	0.994	0.934	0.327
Fold 3 (test)	0.771	0.572	0.758	0.767	0.593	0.802
Fold 4 (train)	0.996	0.948	0.240	0.997	0.951	0.177
Fold 4 (test)	0.618	0.397	0.874	0.642	0.472	0.901
Fold 5 (train)	0.995	0.941	0.231	0.991	0.921	0.380
Fold 5 (test)	0.838	0.699	0.767	0.719	0.524	0.900
Average (train)	0.994	0.934	0.286	0.995	0.941	0.275
Average (test)	0.729	0.540	0.789	0.727	0.532	0.822

Conclusion

SARS-CoV-2 main protease (Mpro) is the most favorable target for COVID-19 drug discovery due to its conservative nature and low similarity with human genes. Structure and binding affinity of protein–drug complexes are of paramount importance for understanding the molecular mechanism in drug discovery. However, there are only two SARS-CoV-2 Mpro inhibitor structures available with binding affinities, highlighting current challenges in COVID-19 drug discovery.

This work presents the reliable binding affinity prediction and ranking of 137 Mpro–inhibitor crystal structures that have no reported experimental binding affinity. We first curate a set of more than 600 Mpro inhibitors with binding affinities from public resources, such as PDBbind, ChEMBL and the scattered literature. Among these inhibitors, 141 are retained based on their high similarity with available Mpro–inhibitor complex structures and built with three dimensional (3D) poses using our MathPose.[6] Together with 44 another SARS-CoV or SARS-CoV-2 Mpro–inhibitor complexes, we compose a training set of 185 reliable SARS-CoV-2 Mpro–inhibitor complexes. Our earlier MathDL models are reformulated with algebraic topology to accommodate 119 new complexes and 17 382 complexes from the PDBbind v2019 general set in both single-task and multitask settings, which have never been available before. The resulting MathDL models are rigorously validated via PDBbind v2016 core set benchmark in which it outperforms state-of-the-art models in the literature. Most importantly, our MathDL achieves promising cross-validation accuracies on the SARS-CoV family inhibitors with the averaged Pearson's correlation coefficient as high as 0.73.

Additionally, the present work unveils that Gly143 of Mpro is the most attractive region to form hydrogen bonds, followed by Glu166, Cys145, and His163. There are 71 inhibitors interacting with SARS-CoV-2 Mpro to form covalent complexes. Those covalent bonds are mostly composed between dicarbon monoxide groups in inhibitors and γ-sulfur on Cys145. There are only one non-covalent complex in our top 10 ranked, namely 5rg1. To provide a potential resource for lead molecule design, we employ the BRICS algorithm to decompose all the inhibitors of the prominent binding site on Mpro and obtain 185 unique fragments.

The predicted binding affinities and their ranking of 137 Mpro–inhibitor crystal structures, the bonding analysis, and the fragment decomposition have significantly extended current knowledge and understanding of SARS-CoV-2 Mpro and inhibitor interactions and, thus offered valuable information toward COVID-19 drug discovery.

Conflicts of interest

There are no conflicts to declare.