key: cord-0749321-gaogd253
authors: Tam, Nguyen Minh; Nam, Pham Cam; Quang, Duong Tuan; Tung, Nguyen Thanh; Vu, Van V.; Ngo, Son Tung
title: Binding of inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro
date: 2021-01-13
journal: RSC advances
DOI: 10.1039/d0ra09858b
sha: fc7dbec26b4b3c3bcfd11a3bdc4aaf93282eb16c
doc_id: 749321
cord_uid: gaogd253

SARS-CoV-2 rapidly infects millions of people worldwide since December 2019. There is still no effective treatment for the virus, resulting in the death of more than one million patients. Inhibiting the activity of SARS-CoV-2 main protease (Mpro), 3C-like protease (3CLP), is able to block the viral replication and proliferation. In this context, our study has revealed that in silico screening for inhibitors of SARS-CoV-2 Mpro can be reliably done using the monomeric structure of the Mpro instead of the dimeric one. Docking and fast pulling of ligand (FPL) simulations for both monomeric and dimeric forms correlate well with the corresponding experimental binding affinity data of 24 compounds. The obtained results were also confirmed via binding pose and noncovalent contact analyses. Our study results show that it is possible to speed up computer-aided drug design for SARS-CoV-2 Mpro by focusing on the monomeric form instead of the larger dimeric one.

The novel coronavirus (2019-nCoV or SARS-CoV-2), a member of the Coronaviridae virus family, has been reported to be able to spread among humans. 1 The virus initially appeared in December 2019 in Wuhan, Hubei Province, China. [2] [3] [4] It shares more than 82% identical RNA genome to the SARS-CoV, SARS-CoV-2 severe cases of respiratory syndromes. 5 Although the bat has been thought of as the original reservoir, the intermediate host is still unknown. 6 Moreover, it is known that the SARS-CoV-2 can endure in aerosol for more than 3 hours, 7 which may be a major factor behind the outbreak of COVID-19 pandemic, which has caused several hundred thousands of deaths worldwide. 5 Therefore, the COVID-19 pandemic becomes an urgency for community health, which requires to develop an effective treatment or vaccine immediately.

Coronaviruses genomes occupy ca. 26-32 kb in length which is the largest sequence among RNA viruses. 8, 9 The SARS-CoV-2 genome encodes more than 20 various structural and nonstructural proteins. Particularly, the SARS-CoV-2 main protease (Mpro), 3C-like protease (3CLP), is one of the most important viral enzymes, having more than 96% similarity with SARS-CoV 3CLP. 9,10 SARS-CoV-2 Mpro cleaves nascent polyproteins, which are produced by the translation of the viral RNA. During this process, 11 non-structural polyproteins are auto-cleaved to become polypeptides, which are required for the viral replication and transcription. 9 Therefore, SARS-CoV-2 Mpro turns out to be an attractive target for antiviral drug development since blocking viral protease can inhibit viral replication and proliferation. 10, 11 Numerous investigations following this strategy have been carried out and shown some initial success. [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] However, unfortunately, an effective drug for COVID-19 is still unavailable.

Currently, it should be noted that the time and cost to advance a drug has been signicantly decreased by using the power of computational approaches. [22] [23] [24] [25] Generally, the binding free energy, DG, between a ligand and an enzyme can be probed via computational approaches. The DG is associated with the experimental inhibition constant, k i , via formula DG bind ¼ RT ln(k i ), where R is gas constant, T is absolute temperature, and k i is a critical metric revealing the nature of binding between two biomolecules. 22 Accurate assessment of the ligandbinding free energy is very important in computer-aided drug design (CADD) problem. 26 In addition, it should be noted that the dimer was shown to be the biologically active form of the SARS-CoV Mpro instead of the monomeric one. 27 Moreover, the SARS-CoV-2 Mpro possibly acts like the SARS-CoV Mpro due to the dissimilitude of only is 12 of 306 amino acids. However, fortunately, the interface of SARS-CoV-2 Mpro does not contain a ligand-binding pocket, 14 the computational screening potential inhibitors for SARS-CoV-2 Mpro are thus possible to carry out based on the monoclinic structure. 28, 29 However, an important question is raised that what is the difference when we use the monomeric form of SARS-CoV-2 Mpro as inhibitor-screening target instead of the dimeric one to reduce CPU time consumption? In this context, the binding affinity of available inhibitors [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] to the monomeric and dimeric SARS-CoV-2 Mpro was examined via docking and FPL schemes. The affinity of the some inhibitors of SARS-CoV-2 to the Mpro was also evaluated. The high correlation coefficients between computational values of monomeric and dimeric systems suggests that we can use the monomeric form of SARS-CoV-2 Mpro as CADD target instead of the dimeric form. Moreover, the similar of Pearson correlation between computed and experimental metrics of SARS-CoV-2 Mpro monomer and dimer was observed. The obtained results can be benecial to the COVID-19 therapy by speeding up CADD progression.

Three-dimensional structures of the monomeric and dimeric SARS-COV-2 Mpro were downloaded from the Protein Data Bank with ID 6Y2F 14 and 6XBG, 30 respectively. The protonation states of the protease at pH ¼ 7.0 was computed using H++ server. 31 Inhibitor structures were taken from the PubChem database. 32 The ligand protonation state were assessed at pH ¼ 7.0 by using a webserver, http://www.chemicalize.com, which is a tool of ChemAxon. The ligand structure was rst optimized using quantum mechanics (QM) simulation with the B3LYP functional at 6-31G(d) level of basis set.

The binding position and affinity of ligands to the monomeric and dimeric SARS-CoV-2 Mpro were probed via the Autodock Vina package (cf. Fig. 1 ). 33 In particular, the AutodockTool 1.5.6 was manipulated to topologize the ligands and receptors. 34 The docking parameter was selected referring to the previous study, [35] [36] [37] in which the exhaustiveness is of 8. The obtaineddocking result was chosen as the highest binding affinity conformations. The grid center was selected as the geometric center of the a-ketoamide 13b and UAW246 compounds, which correspond to the monomeric and dimeric Mpro, respectively. 14, 30 The grid size was chosen as 24 Â 24 Â 24Å, which entirely cover the ligand-binding cle of the Mpro. 36, 37 Steered-molecular dynamics simulations GROMACS version 5.1.3 (ref. 38 ) was used to simulate the solvated complex involving the ligand and monomeric/dimeric SARS-COV-2 Mpro. The protease and inhibitor were topologized via the Amber99SB-ILDN 39 and general Amber force eld (GAFF), 40 correspondingly. It should be recorded that the ligand parameterization was completed by using AmberTools18 and ACPYPE approaches. 41, 42 In particular, the inhibitor atomic charges were assigned via the Restrained Electrostatic Potential (RESP) method 40 through QM investigation at the level of B3LYP/6-31G(d,p). During QM simulation, the implicit water model, 3 ¼ 78.4, was involved. The monomeric and dimeric SARS-CoV-2 Mpro + inhibitor were inserted into a rectangular periodic boundary condition (rPBC) box with a dimension of (9.8, 5.9, 8.7) and (9.4, 9.0, 12.1) nm, respectively. The corresponding box volumes of the monomeric and dimeric systems are 506.28 and 1013.82 nm 3 , respectively. Therefore, the total atoms of these systems approximately are 50 000 and 100 000 atoms, respectively.

The atomistic simulation was performed utilizing the parameters referred to the prior works. 36, 37 Particularly, the MD time step is 2 fs. The noncovalent pair was affected within a radius of 0.9 nm. The electrostatics interaction was assessed implementing the fast particle-mesh Ewald electrostatics scheme. 43 The SARS-CoV-2 + inhibitor was then optimized and equalized throughout the EM, NVT, and NPT imitations. The NVT and NPT imitations were operated during intervals of 0.1 and 2.0 ns, correspondingly. Moreover, the SARS-CoV-2 Mpro Ca atoms were restrained during these imitations via a small harmonic force with a value of 1000 kJ mol À1 nm À2 per proportions. The relaxed conformation of the SARS-CoV-2 Mpro + inhibitor was then employed as initial structure of FPL simulation. During which, the inhibitor was pulled out of the binding cle under effect of an externally harmonic force with parameters of k ¼ 0.005 nm ps À1 and v ¼ 600 kJ mol À1 nm À2 for pulling speed and cantilever spring constant (cf. Fig. 1 ), respectively. 37 Totally, 8 independent trajectories were carried out to assess the ligand-binding affinity. 

An intermolecular nonbonded contact was enumerated when the minimum distance of nonhydrogen atoms of a residue to the inhibitor was smaller than 0.45 nm. A hydrogen bond was enumerated when the angle between donor, D, hydrogen, H, acceptor, A, is larger than 135 and the distance between D and A is smaller than 0.35 nm. The error of computations was computed through 1000 rounds of the bootstrapping method. 44

The preliminary investigations of binding pose and affinity of the trial inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro were initially estimated by a molecular docking method. Autodock Vina, 33 a very efficient molecular docking approach with a successful-docking rate up to 81%, 35 would be able to complete this task. We have thus docked 24 available inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro using Autodock Vina referring to the previous study. 36, 37 Moreover, the successful-docking rate of Autodock Vina on the SARS-CoV-2 Mpro is also appropriate due to the recent report with the number is 67%. 45 By using exhaustiveness 8 as suggested in the previous work, 35 the results were rapidly obtained in few hours (Tables 1 and S1 of the ESI le †). Interestingly, in good agreement with the previous studies, 45 Autodock Vina adopts an appropriate results compared with the experimental affinities, in which, the correlation coefficient of the monomeric target, R Monomer Dock ¼ 0.50 AE 0.15, is similar to that of the dimeric target, R Dimer Dock ¼ 0.50 AE 0.12. It should be noted that the obtained correlation is in good agreement with the previous benchmark of Autodock Vina over 800 various complexes, R exhaustiveness 35 Moreover, the root mean square error (RMSE) between calculated and experimental values are small indicating the accurate result, in which, the monomeric system gives a value of RMSE Monomer Dock ¼ 1.89 AE 0.15 kcal mol À1 and the dimeric system adopts a metric of RMSE Dimer Dock ¼ 1.65 AE 0.17 kcal mol À1 . Furthermore, comparing the docking results of ligands to monomeric and dimeric Mpro would reveal the correlation of using monomeric and dimeric systems as docking target. The obtained results indicate that the docking of ligands to monomeric Mpro well correlates with the corresponding values of the dimeric system, R Monomer-Dimer Dock ¼ 0.85 AE 0.08 (Fig. 2) . The accuracy of docking results in different targets is also high due to the small value of the obtained RMSE, RMSE Monomer-Dimer Dock ¼ 0.66 AE 0.10 kcal mol À1 . It may thus argue that we are able to use the monomeric form of SARS-CoV-2 Mpro instead of the dimeric form for screening potential inhibitors of the Mpro via the docking approach. It should be noted that the computed error bars was obtained via 1000 rounds of the bootstrapping method. 44 Furthermore, the binding pose of inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro is in good agreement together since espousing the root-mean-square deviation (RMSD) of 0.21 AE 0.02 nm (cf. Fig. 3 and Table S1 of the ESI †). It should be noted that the RMSD of the ligand-binding poses, which is smaller than 0.20 nm, normally counted as the conformations locating in the same cluster. 35, 46 The difference in docking poses probably causes the uncorrelation between docking energies of monomeric and dimeric systems. The structural observation is thus conrmed the obtained docking energy above. The molecular docking with larger exhaustiveness, which selected as 56 and 400 according to the previous study, 35 were also performed in order to validate the convergence of the docking scheme. In total we used three different values of exhaustiveness including 400, 56, and 8 which are denoted as long, medium, and short options, respectively. The associations of the docking simulations for monomer and dimer with respect to experiment are shown in Table 1 . Interestingly, changing the docking exhaustiveness parameter from short to medium and/or long does not have a signicant impact on the correlation coefficient and RMSE, which is consistent with the prior benchmark. 35 (Fig. 4) , respectively. Overall, the docking simulations provide similar results when ligands docked to monomeric and dimeric systems and the default option of Autodock Vina is appropriate for performing docking simulations.

As mentioned above, the binding affinity of 24 available inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro was appropriately probed using molecular docking calculations. However, it should be noted that the dynamics of receptors were not considered in docking simulations, and the number of trial docking poses was restricted. To overcome this limitation we have operated the atomistic simulations which serve as a validation for the docking results. [47] [48] [49] Moreover, FPL is an efficient computational approach to assess ligand-binding affinity with a suitable time-consuming calculation. 50,51 Furthermore, the scheme was successfully applied to the monomeric SARS-CoV-2 Mpro system recently. 36, 37 The FPL approach is thus used to probe the binding affinity of 24 available inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro. In the simulations, the ligand binding pose was optimized over short canonical and isothermal-isobaric simulations. The equilibrated ligand was then pulled to translocate from bound to unbound states. The maximum of pulling force, called rupture force, and pulling work are typically assumed to correlate with ligand-binding affinity. It should be noted that the rupture force corresponds to the point that the non-covalent bond contact between a ligand and a receptor was terminated.

The computed values of the rupture force and pulling work were shown in Table 2 . The denoted pulling force and work proles were described in Tables S2 and S3 of the ESI le. † The shape of both pulling force and work appear reliable when compared to the previous exertion. 50, 51 In particular, starting at zero, the pulling force quickly increases to the maximum value, then suddenly drops to zero due to the loss the noncovalent bond contact to the receptor. During this process, recorded-pulling work speedily rises from zero value to a stable value, corresponding to the distance at which the contact between protein and inhibitor is vanished. Moreover, the rupture force F Monomer Max of monomeric Mpros diffuses in the range from 295.0 to 977.6 pN corresponding with the spreading of pulling work W Monomer from 13.7 to 106.1 kcal mol À1 . Besides that, the matching metrics of dimeric Mpros forms in the range from 307.7 to 860.2 pN and 22.4 to 85.7 kcal mol À1 , correspondingly. It should be noted that the computed works are signicantly larger than the magnitude of experimental binding affinity, which diffuses in the range from 6.39 to 10.13 kcal mol À1 , since applied large cantilever and high pulling velocity. 50 Although the discrepancy can be reduced to zero by using a small cantilever and an extremely low pulling velocity, it is not appropriate since it requires to perform several trajectories with hundred nanoseconds each. 52 Furthermore, previous investigations revealed that although reducing the magnitude of cantilever spring constant and pulling velocity was able to enlarge the accuracy of the estimations, the observed results are approximately the equivalent as those at high pulling velocity. 50 a The experimental binding free energies were gained based on IC 50 value, [12] [13] [14] [15] [16] [17] [18] approximating that the one equals to the inhibition constant k i . The unit of force and energy/work are in pN and kcal mol À1 , respectively.

In practice, the rupture force has been used as a predictor of ligand-binding affinity based on the assumption that a ligand binds with a higher affinity requires a stronger pulling force to dissociate it from binding cle. 53 Using the rupture force as a proxy to ligand-binding affinity, numerous investigations were successful in predicting the ligand-binding affinity to various targets. 53, 54 Here, the average of rupture forces were estimated over 8 independent FPL trajectories (cf. Table 2 ). The correlation coefficient to experiments, obtained results of monomeric systems, is R Monomer Force ¼ À0.69 AE 0.10; while the analogous value of dimeric forms is R Dimer Force ¼ À0.60 AE 0.15. Because the correlation coefficients appear to be the same within the error range, we may conclude that there is no difference when using monomer or dimer as a CADD target. Moreover, the correlation between rupture forces obtained over SARS-CoV-2 Mpro monomer and dimer is appropriate with a value of R Monomer-Dimer Force ¼ 0.67 AE 0.10 (cf. Fig. 5 ). It may be argued that the recorded rupture forces of the SARS-CoV-2 Mpro monomer is quite similar to the dimeric one, although the correlation is smaller than obtained values via docking approach.

The work of pulling force was assessed via formula

where v is pulling velocity and F(t) is pulling force. In isothermal-isobaric simulations, W is related to the experimental binding affinity via Jarzynski equality. 55 Therefore, utilizing W to estimate the ligand-binding affinity commonly acquires a better accurate result in comparison to rupture force. 47, 50, 54 The obtained results reaffirmed this statement. The correlation coefficients of the monomeric and dimeric SARS-CoV-2 Mpro are R Monomer Work ¼ À0.69 AE 0.09 and R Dimer Work ¼ À0.65 AE 0.13, respectively. Although, the computational accuracy targeting the SARS-CoV-2 Mpro monomer is slightly larger than that of the dimeric system, the difference in correlation coefficients is insignicant implying that the monomeric form of SARS-CoV-2 Mpro can be used as CADD target instead of the dimeric one. Moreover, it should be noted that FPL outcomes based on the short NPT simulation of only 2.0 ns. The results possibly limited since the complex may not gain the equilibrium states as reported in the recent work. 45 The association of computed pulling works of the monomeric and dimeric SARS-CoV-2 Mpro was probed and shown in Fig. 6 . Over the bootstrapping examination, the correlation coefficient is R Monomer-Dimer Work ¼ 0.78 AE 0.06 conrming the observation above. The non-enhancement of FPL outcomes compared with docking results possibly occurs due to the increase of the different binding poses to various targets, which was 20, 21 approximating that the one equals to the inhibition constant k i . The unit of force and energy/work are in pN and kcal mol À1 , respectively. described below. Overall, it may be concluded that we can perform the inhibitor screening for SARS-CoV-2 Mpro with smaller computing resources since targeting the monomeric form.

In addition, the MD-rened ligand-binding affinity results are conrmed since the RMSD between ligand-binding poses to the monomeric and dimeric forms is of 0.23 AE 0.02 nm only. It should be noted that the RMSD metrics were calculated based on the last snapshot of NPT simulations, which were utilized for the binding free energy prediction via the FPL scheme. Moreover, the RMSD of MD-rened structure was slightly larger than docking results due to the effects of the conformational entropy. Moreover, the intermolecular hydrogen bond analyses suggests that three residues including Thr25, Asn142, Gly143, and Glu166 are critical residues controlling the binding mechanism of the inhibitors to both monomeric and dimeric SARS-CoV-2 Mpro (cf. Table S4 of the ESI †).

The binding affinity of 19 available SARS-CoV-2 inhibitors to the protease was also computed. The obtained results were reported in Tables 3 and S5-S7 of the ESI le. † Although IC 50 of these inhibitors was taken from cell-culture experiments 21 implying that drug targets possibly differ from the protease as well as RNA polymerase, appropriate correlations between computed outcomes and respective experiments were obtained. Therefore, it may be concluded that many compounds inhibit SARS-CoV-2 Mpro. In particular, Autodock Vina adopts a larger correlation coefficient on the monomer, R Monomer Vina z 0.84, than dimer, R Dimer Vina z 0.52, whereas the R Monomer-Dimer Dock z 0.70 implies the consistency of our observation. The SARS-CoV-2 Mpro + inhibitor structure provided by the docking approach was then used as the starting shape for FPL calculation. The association between computed values of the monomeric and dimeric Mpro were obtained. The disparity attained due to some compounds possibly aiming on the polymerase rather than the protease. 19 Conclusions Both of Autodock Vina and FPL simulations were conrmed to be able to appropriately estimate the ligand-binding affinity of the SARS-CoV-2 Mpro in both monomeric and dimeric forms. The assessed results suggested that the monomeric form of SARS-CoV-2 Mpro can be used as a CADD target instead of the dimeric form. In particular, the high correlation coefficients between computational ligand-binding affinities of the monomeric and dimeric SARS-CoV-2 Mpro via docking and FPL simulations are obtained, which are R Monomer-Dimer Dock ¼ 0.85 AE 0.09 and R Monomer-Dimer Work ¼ 0.78 AE 0.06, respectively. Moreover, the correlation coefficient between the rupture forces of two targets is roughly appropriate with a value of R Monomer-Dimer Force ¼ 067 AE 0.10. It should be noted that the observation is in good agreement with structure analyses with the RMSD between ligand-binding pose to the monomeric and dimeric forms of 0.21 AE 0.02 and 0.23 AE 0.02 nm for docking and MD-rened structures, respectively. Furthermore, the appearance of different binding poses possibly results in the non-perfect correlation between the computational values of two targets.

The obtained computational values of the monomeric and dimeric SARS-CoV-2 Mpro correlate to experiments in similar amounts. The docking approach formed a Pearson correlation of ca. R Dock ¼ 0.50 to both targets. FPL approach enhanced the accuracy of the calculated ligand-binding affinity since a correlation is ca. R Work z À0.65 over both receptors. The accuracy of FPL simulations probably increases due to extending the NPT simulation time as reported in the recent work. 45 In addition, in good agreement with the previous observation, 35 the molecular docking by Vina package rapidly converged since the correlation coefficient between computed and experimental values did not change when the docking option was altered. The RMSE of docking results also unchanged upon these alterations. Finally, it may be concluded that for SARS-CoV-2 Mpro system the pulling work is better than rupture force in predicting the ligand-binding affinity. It is well compatible with earlier probe various protein-ligand complexes. 47, 50, 54 

There are no conicts to declare.

Potential Non-Covalent SARS-CoV-2 3C-like Protease Inhibitors Designed Using Generative Deep Learning Approaches and Reviewed by Human Medicinal Chemist in Virtual Reality

Antibiotics: Methods and Protocols