key: cord-0755113-vxeuqzdm authors: Ma, Yufei; Tao, Yulian; Qu, Hanyang; Wang, Cuihong; Yan, Fei; Gao, Xiujun; Zhang, Meiling title: Exploration of plant-derived natural polyphenols toward COVID-19 main protease inhibitors: DFT, molecular docking approach, and molecular dynamics simulations date: 2022-02-14 journal: RSC advances DOI: 10.1039/d1ra07364h sha: f95c98e1c2ba32ef25ac6a6584b8383d35805eac doc_id: 755113 cord_uid: vxeuqzdm Recent outbreaks of coronavirus have brought serious challenges to public health around the world, and it is essential to find effective treatments. In this study, the 3C-like proteinase (3CLpro) of SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) has been considered as an important drug target because of its role in viral replication. We initially optimized 251 compounds at the PM7 level of theory for docking with 3CLpro, and then we selected the top 12 compounds for further optimization with the B3LYP-D3/6-311G** method and obtained the top four compounds by further molecular docking. Quantum chemistry calculations were performed to predict molecular properties, such as the electrostatic potential and some CDFT descriptors. We also performed molecular dynamics simulations and free energy calculations to determine the relative stability of the selected four potential compounds. We have identified key residues controlling the 3CLpro/ligand binding from per-residue based decomposition of the binding free energy. Convincingly, the comprehensive results support the conclusion that the compounds have the potential to become a candidate for anti-coronavirus treatment. The COVID-19 acute respiratory tract disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) was rst reported in December 2019 in Wuhan, China. [1] [2] [3] The disease was declared a pandemic by the World Health Organization (WHO) on March 11, 2020. 4 In the absence of potential and specic therapeutics, COVID-19 cases were rapidly increasing, and over 232 million cases with more than 4.75 million deaths were reported up to September 2021. 5 To date, though the development of such drugs is underway, 6, 7 there are still no commonly used effective drugs for treating such CoVs. SARS-CoV-2 is a member of the Coronaviridae family, the same family for other beta coronaviruses, such as Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) and Middle East Respiratory Syndrome Coronavirus (MERS-CoV). [8] [9] [10] Recent studies have shown that SARS-CoV-2 shares $89% sequence similarity with other SARS-CoVs. 8 SARS-CoV-2 contains two kinds of proteins, i.e., 4 structural proteins (spike, envelope, membrane, and nucleocapsid) and 16 nonstructural proteins. In the autoproteolytic cleavage catalyzed by PLpro enzymes (papain-like cysteine protease) and 3CLpro enzymes (3C-like proteinase), SARS-CoV-2 genetic material enters the host cell and borrows ribosomes to translate into 16 nonstructural proteins. Because of its mechanistic signicance, 3CLpro, also known as Mpro, is a central target for the effective inhibitors (antiviral drugs) against SARS-CoV-2 and other known CoVs. 11 The crystal structure of main protease 3CLpro of SARS-CoV-2 in the complex with an inhibitor N3 (PDB code 7BQY, shown in Fig. 1 with PyMOL 12 ) has been revealed, with the resolution of 1.7Å. 13 This protease has been considered to be a promising therapeutic target for COVID-19 since 3CLpro plays an important role in the viral replication process while no close homolog exists in humans. 10, 14 Natural compounds have been used for medical treatments since ancient times because of their low cost, minimum sideeffects, and abundant therapeutic resources. In some previous studies, the compounds based on natural products were discovered as inhibitors against viruses, such as polio-virus type 1, the RNA-dependent RNA polymerase (RdRp), respiratory synchronization (sin-SISH-Uhl), and most notably as purinebased inhibitors of viral replication complex component. 15, 16 The avonoids are phenolic compounds that are among the most plentiful and prevalent natural compounds. There are ten major sub-groups of avonoids, i.e., aurones, biavonoids, catechins, avanones, avanols, chalcones, avanonols, avans, avones and isoavones. 17 Baicalin (a known avonoid) has an EC50 value of 12.5-25 mg ml À1 at 48 h without causing signicant cytotoxicity. 18 Flavonoids luteolin and quercetin have the ability to prevent SARS-CoV from infecting host cells. 19 Gallocatechin gallate (GCG) display good inhibition toward 3CLpro with IC50 values of 47 mM. 20 Because of these properties, avonoids could be potential candidates to interfere with the life cycle of coronaviruses. 21, 22 According to in vtiro studies, epigallocatechin gallate (EGCG) inhibited 3CLpro by 85% at a concentration of 200 mM, and had an IC50 measure of 73 AE 2 mM. 20 The anthocyanin of red-eshed potato inactivated both inuenza viruses A (IVA) and B (IVB). 23 Hence, avonoids, as other natural products, continue to be a promising source of anti-coronavirus therapeutics. Traditional Chinese medicine has a long history and has accumulated rich experience in the use of natural medicines. 24 Aer the careful literature search on the topics related to avonoids and viral diseases (specically the coronavirus diseases), we established our screening library with a total of 251 kinds of polyphenols included. In light of the long and expensive process of developing compounds into approved drugs and the rapidly evolving nature of the epidemic, a combination of computational methodologies is a promising approach for identifying potential drug candidates from compound libraries. 25, 26 In the present paper (as shown in Fig. 2) , rstly, we optimized the molecules in our screening library. In order to further study the composition characteristics of the molecules, the compounds based on docking scores were studied by the quantum mechanics method. We also used molecular dynamics simulation to identify the stability of promising compounds with Mpro. Then, we performed binding energy calculations using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method to evaluate the binding affinity of the compounds and key residues. In silico studies were performed for identifying bioactive compounds that can inhibit COVID-19 Mpro effectively, as well as providing useful information for future studies. The X-ray crystal structure of 3CLpro in complex with the inhibitor N3 (PDB ID 7BQY) with the resolution of 1.7Å was used for our study, and the crystal structure was prepared by removing the N3 ligand and adding hydrogen atoms, as well as computing the Gasteiger charge using the AutoDock v4.2 program. 27 Our research focuses on avonoids and viral diseases (specically Coronavirus diseases), and in order to broaden the scope of our search we collected compounds related to avonoids from Chinese herbs. The SDF structures of the selected 251 ligands (shown in Table S1 †) were retrieved from the Pub-Chem database (https://pub-chem.ncbi.nlm.nih.gov/). We adopted the PM7 method implemented in Gaussian 16 (ref. 28) to optimize these 251 compounds, and the optimized geometries were used to do docking study with the Autodock Vina program. 29 Based on the scores of docking, we screened 12 compounds out of the 251 candidates for the further geometry optimization with the B3LYP-D3/6-311G** method, which is more reliable than PM7. Finally, a second docking for the 12 compounds was performed with Autodock Vina, and 4 top-score compounds were selected. The ligands were saved and converted to PDBQT format using Raccoon. 30 The top 4 of 12 highest binding energy and best-docked conformation were considered for the further quantum chemistry calculations. A grid box with dimensions 60 Â 60 Â 60 centered at N3 in protein (X ¼ 9.269, Y ¼ À1.127, and Z ¼ 22.118) was used as the search area. We used the soware Discovery Studio Visualizer 31 to illustrate the interactions between the best four ligands and proteins, including hydrogen bonds, atom accessibility, and hydrophobic interactions (Fig. 3 ). In order to predict molecular properties, we have carried out the related quantitative calculation. 32 Among various computational methodologies, DFT has become the preferred method for investigating molecular properties, which has been extensively used to compute the electronic structure properties. Because of comparatively better results and less computational cost, 33 we choose the DFT/B3LYP methods 34 for the optimization of the screened compounds against SARS-CoV-2 Mpro using Gaussian 16 soware. Electrostatic Potential (ESP) 35 were generated for the selected molecules with Multiwfn, 36 and the ESP surfaces are shown by VMD 37 to analyze the electrostatic induced weak interaction. According to the theory of the Conceptual Density Functional Theory (CDFT), we calculated the chemical parameters to show more molecular properties within the framework of the HSAB theory, including E HOMO , E LUMO , ionization potential (IP), electron affinity (E A ), energy gap (E g ), electronegativity (c), chemical potential (m), chemical hardness (h), chemical soness (S), electrophilicity (u) and nucleophilicity (N). The ionization energy and electron affinity are related to the HOMO and LUMO energies, respectively, based on Koopmans' theorem. 38 The related formulas are as follows: The four ligands with the best docking score were selected for the MD simulation using Gromacs 39 to evaluate their intermolecular interactions with 3CLpro and the stability of the complexes. The ligands were parameterized with the general AMBER force eld (GAFF) 40 using the Acpype program. 41 The restrained electrostatic potential (RESP) approach at the B3LYP/ 6-311G** level was applied to calculate the atomic partial charges of the analogs with the assistance of Gaussian16 soware. 42 Protein topology and coordinate les were generated using the Amber99SB force eld 43 provided in Gromacs. 39 The protein-ligand complex was contained in a dodecahedron and solvated with an explicit TIP3P 44 water model within a periodic boundary box of distance 1.0 nm xed in between the protein and dodecahedron box, and the Na + ions were added to neutralize the solvated system that followed by quick energy minimization with the steepest descent minimization algorithm. This was followed by a restrained constant number of particles, volume, and temperature (NVT) ensemble equilibration for 500 ps and then simulated for 1 ns in the NPT ensemble at 310 K (ref. 45 ) and 1.0 bar. 46 The particle mesh Ewald method was used to calculate the long-range electrostatics. 47 Modied Berendsen thermostat and Parrinello-Rahman barostat were used for temperature and pressure coupling, respectively. Finally, unrestrained 100 ns production simulations were carried out for the systems at 310 K and 1 bar atmospheric pressure. For MD simulation analysis, the results of RMSD and bonding interactions were performed using Gromacs tools. 48 Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) has been regarded as a competitive method in the eld of binding free energy calculation, 49 and is applied to estimate the relative binding free energies for the four nal compounds and to evaluate their relative stability of the biomolecular structures against SARS-CoV-2 in our studying. 50, 51 We use g_mmpbsa that was developed to enable the use of the MM-PBSA method in conjunction with the Gromacs package to do the calculation. 52 In general terms, the binding free energy of the protein with ligand in solvent can be expressed as: where, G complex is the total free energy of the protein-ligand complex, and G protein and G ligand are total free energies of the isolated protein and ligand in solvent, respectively. g_mmpbsa can also be used to estimate the energy contribution per residue to the binding energy. To decompose the binding energy, at rst DE MM , DG polar and DG non-polar were separately calculated for each residue and were then summed up to obtain the contribution of each residue to the binding energy. Molecular docking is an effective tool in drug design, 53 which speeds up this process. We use molecular docking technology to screen compounds based on their docking scores and identify the key interactions of ligands within the active site of a target protein. 53 The docking scores of the rst 12 compounds were shown in Table 1 . Based on docking energy values of À9.7, À9.7, À9.5, and À9.5 kcal mol À1 , respectively, the four compounds (isocryptomerin, hinokiavone, amentoavone, and neocryptomerin) with favorable binding affinities in the active pocket of SARS-CoV-2 were selected for the further studying. The observed hydrogen bond interactions and hydrophobic interactions were shown in Fig. 3 . Neocryptomerin ts well the active cavity of 3CLpro and is stabilized by three hydrogen bonds (at the sites of Cys145, Ser46, and Thr24). Besides, some residues like Met165, (Pi-alkyl interaction), Met49 (Pi-sigma interaction), Cys145 (Pi-sulfur interactions), His41 (Pi-Pi Tshaped interactions), and Thr25 (Pi-donor hydrogen bond interactions) are found to undergo hydrophobic interactions. Amentoavone forms Hydrogen bonds with Asn142, Ser46, and Thr24, while forming multiple interactions with His41 and Met165. Four hydrogen bonds of 3CLpro/hinokiavone are formed within the active site, while His41 (Pi-Pi T-shaped interactions), Cys145 (Pi-sulfur interactions), and Met165 (Pialkyl interactions) also interact with the compound. The main residues of isocryptomerin forming hydrogen bonds with proteins include Cys145, Ser46, Thr45, and Thr24, while other interactions are also spotted, such as Pi-sulfur, alkyl, Pi-Pi Tshaped, Pi-donor hydrogen bond, and Pi-lone pair interactions between the isocryptomerin and Mpro. Aer the docking screening, four compounds were found to have substantial binding energy toward the 3CLpro, and are expected to inhibit the 3CLpro activity, thus being the candidates for further study. Studying the ESP of various compounds is useful for nding new chemical or biochemical compounds that can be used in the discovery and design of new therapeutics and for developing them into useful medicines. The interactions between biomolecules and compounds mainly include hydrogen bonds and halogen bonds dominated by electrostatic interactions, and ESP is the real space function that is most closely related to the electrostatic effect, and it is very suitable to analyze the electrostatic-induced weak interaction. Therefore, the analysis of electrostatic potential is helpful to understand the chargedependent properties in the molecular structure of the complex. 35, 54 Fig. 4 shows the ESP distribution of the molecule as well as the active sites and comparative reactivity of the atoms. The negative region is presented in blue, and the positive region in red. The histogram in Fig. 4 shows that the large portion vdW surface of them have ESP value within À10 to +10 kcal mol À1 . Only a tiny part of the vdW surface has ESP value larger than +60 kcal mol À1 or less than À65 kcal mol À1 . Even small areas with minima and maxima of ESP can be very important to establish hydrogen bonds or electrostatic interactions. The most negative region corresponds to the oxygen in 4-Oxo-4h-1-benzopyran, and the most positive region to the hydrogen in hydroxyl. Surface local minima and maxima of ESP for neocryptomerin, isocryptomerin, hinkiavone and amen-toavone are labeled by green and orange spheres, respectively, with their values shown in Fig. 4 . These extremes are the nucleophilic and electrophilic sites, which are most likely to be active with residues. As expected, the ESP results are consistent with the hydrogen bond and polar interactions depicted in Fig. 3 . According to the frontier molecular orbital theory, the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO) are the keys to determining the chemical reaction of the system. 55 Therefore, the study of frontier orbitals of the drug molecules can provide important information for the exploration of the mechanism of action and the identication of active sites. 56 HOMO-LUMO energy gaps of similar systems indicate the stability of molecules, with a smaller energy gap representing a more unstable molecule. 57 The HOMO and LUMO orbitals for the screened compounds and their HOMO-LUMO energy gaps are shown in Table 2 and Fig. 5 , and LUMO-HOMO energy gaps decreased in the following order: amentoavone < isocryptomerin < hinokifavone < neocryptomerin. The smaller the gap is, the easier it interacts with the protein. As we all know, the quantum chemical parameters of the compound are very important for their practical applications. Therefore, for the selected four potential inhibitors, we calculate some physical and chemical properties dened by the Conceptual Density Functional Theory (CDFT), 59 such as E HOMO , E LUMO , ionization potential (IP), electron affinity (E A ), electronegativity (c), chemical potential (m), chemical hardness (h), chemical soness (S), electrophilicity (u) and nucleophilicity (N). Compounds can be studied using CDFT to predict reactive sites and the reactivity, as well as their properties. 38, 60 By comparing the parameters we calculated, we can further screen the promising compounds. Another functional (M06-2X) has also been included for comparison (shown in Table S2 in ESI †), and the CDFT results are comparable with those from the B3LYP. Amentoavone is more reactive than others because of its low chemical hardness (h) and naturally high soness (S) based on Pearson's HSAB principle. 61 From Table 2 , since the hardness values h of amentoavone (6.42 eV) and hinokiavone (6.46 eV) are smaller than those of neocryptomerin (6.47 eV) and isocryptomerin (6.48 eV), amentoavone and hinokiavone are less stable or more reactive than two others. Electronegativity (c) and chemical potential (m) measure the ability of an atom or a functional group to attract electrons. 62 These two values in Table 2 indicate that amentoavone has the most electron attracting ability, while hinokiavone possesses the least electronegativity. The electrophilicity index denes the tendency of an electrophile to acquire a given amount of electron density and the resistance for a molecule to exchange electron density with the environment. 63 Thus, among all compounds, amento-avone shows higher electrophilicity at 1.27 eV than other three compounds, and hinokiavone shows only less value at 1.21 eV. Thus, amentoavone and hinokiavone are more active than neocryptomerin and isocryptomerin in terms of polar reactions activity. The screened compounds are characterized using quantitative chemical parameters in this part. Calculations were done to establish relevant data so that we can gain a better understanding of the properties of the compound. In order to evaluate the stability of the initial structure of the protein and ligand complex, we proceed with the molecular dynamics simulation. In the MD simulation, the root-meansquare deviations (RMSDs) of backbone atoms relative to their respective initial positions are observed (Fig. 6) , and the average RMSD values for neocryptomerin, isocryptomerin, hinoki-avone, and amentoavone are 0.16 nm, 0.15 nm, 0.23 nm, and 0.14 nm, respectively. The highest RMSD deviation (0.33 nm) is obtained for the 3CLpro/hinokiavone system, while the least RMSD (0.09 nm) is obtained for 3CLpro/amentoavone. As can be seen in Fig. 6 , all compounds tend to the equilibrium state aer 80 ns. Neocryptomerin, isocryptomerin, and amentoavone display an average RMSD of <0.2 nm. However, hinokiavone shows higher uctuations as compared to others, and an average RMSD of <0.3 nm was noted from Fig. 6 . 3CLpro/ neocryptomerin has a slight acceptable uctuation at 50 ns, and this system is stable throughout the simulation. For 3CLpro/isocryptomerin, a slight disturbance between 50 and 55 ns is shown, but it remains stable during the simulation. 3CLpro/hinokiavone, which uctuates slightly before 60 ns, gains stability and enters the production stage aer 60 ns. In the case of 3CLpro/amentoavone, it uctuates slightly at about 75 ns, but it is stable at other time periods. The data in Fig. 6 indicate that all these screened compounds are stable in the entire simulation process except for some negligible small uctuations, proving that these compounds are able to bind steadily to the proteins without affecting the overall topology of SARS-CoV-2 Mpro. The binding free energies and energy decomposition analyses for the four screened compounds are evaluated with the help of Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach. 50 Through the production process, the coordinates and energy values are gathered every 100 ps. It provides individual components, such as van der Waals energy, electrostatic energy, polar solvation energy, non-polar solvation energy, and total binding free energy of all the systems. In Fig. 7 , the components of the binding free energy for the inhibitors in the presence of 3CLpro are plotted, and the corresponding data are presented in Table 3 . This study indicates that the van der Waals attraction, electrostatic interactions (except for isocryptomerin), and non-polar solvation energy are the main contributions to the binding, whereas the polar solvation free energy weakens the complexation. Table 3 shows that the estimated binding free energies for the 3CLpro/amentoavone, 3CLpro/hinokiavone and 3CLpro/neocryptomerin complexes are À20.56, À20.50, and À13.10 kcal molÀ1, respectively, while that of 3CLpro/ isocryptomerin complex is positive (+0.30 kcal molÀ1). In short, the binding affinity to 3CLpro decreases in the following order: amentoavone > hinokiavone > neocryptomerin > isocryptomerin. It is further revealed in Table 3 that the van der Waals attraction and non-polar solvation energy of the 3CLpro/ neocryptomerin complex are almost the same as those of the 3CLpro/amentoavone complex and the 3CLpro/hinokiavone complex, but its electrostatic energy (+0.35 kcal mol À1 ) is positive, different from other complexes. Thus, the total binding energy is less. For the 3CLpro/isocryptomerin complex, the van der Waals energy (À7.66 kcal mol À1 ), which plays an important role in stability, allows its binding energy to be positive. Amentoavone has the highest binding energy and is relatively stable, while hinokiavone has its binding energy very close to amentoavone. From Table 3 , it can be seen that the stability- promoting components for the two compounds (amentoavone and hinokiavone) are van der Waals energy, electrostatic interaction, and non-polar solvation energy, while the stabilityinhibiting component is polar solvation energy. The binding energies for other complexes are also calculated and shown in Table S3 in ESI. † Overall, amentoavone, hinokiavone, and neocryptomerin can be considered as the leading compounds of rational compounds against COVID-19. In order to identify the crucial residues responsible for stabilizing the interactions between the inhibitors and 3CLpro, binding free energy contributions from individual residues were decomposed with the MM-PBSA method. Decomposing the binding free energy on a per-residue basis into contributions from internal energy, polar solvation energy, and nonpolar solvation energy for residues, as shown in Fig. 7 . Key residues with a cut-off value of À0.5 kcal mol À1 were identied from the free energy binding spectra. It is evident from Fig. 8 In the case of the 3CLpro/isocryptomerin, there are two residues VAL-125 (À0.28 kcal mol À1 ) and TYR-126 (À0.29 kcal mol À1 ) with binding energy values higher than À0.25 kcal mol À1 . This might be the reason for the poor affinity of isocryptomerin against 3CLpro. The results showed that MET-49 and MET-165 are the key residues, which provided a clue for further research. Overall, the identication of hotspot residues from our analysis can facilitate the discovery of new selective inhibitor against COVID-19. Computational medicinal chemistry has become an important element of modern drug research. The sudden outbreak of the new coronavirus has caused panic all over the world and seriously threatened the health of human beings. However, there is no effective drug now, and the discovery of an effective drug is imminent. Computer-aided drug discovery (CADD) methodologies have emerged as powerful tools in the drug discovery process and could reduce the cost of research and, most importantly, speed up drug discovery. In order to have a more comprehensive understanding of the properties of the drug candidates, we not only did molecular dynamics research but also calculated their molecular properties through quantum chemistry calculations. First of all, we screened 12 candidates from 251 compounds through molecular docking, and the rst 4 compounds were selected aer further structural optimization. These selected four molecules were calculated by the DFT/ B3LYP methods to predict LUMO/HOMO gap, ESP, and other properties related to molecular interacting abilities (E HOMO , E LUMO , ionization potential, electron affinity, electronegativity, chemical potential, chemical hardness, chemical soness, electrophilicity, and nucleophilicity). Aer that, molecular dynamics simulation was carried out to evaluate the stability of the selected compounds, and we have also investigated the hotspot residues controlling the receptor-ligand binding. Based on the results of our study, hinokiavone, and amentoavone are considered as promising compounds against SARS-CoV-2. This is coherent with the previous ndings that these two compounds have antiviral effects. We found that the compound hinokiavone and its derivatives have antagonistic effects against inuenza 64 and dengue virus, 65,66 while amentoavone and its derivatives have inhibitory effects against CVB3, 67 HBA, 68 dengue virus 66 and SARS-CoV. 69 Our ndings suggest that these two compounds might serve as a good starting point for further investigations of potential drugs against SARS-CoV-2. We hope that our ndings in the present paper will be helpful for the experimental study to develop effective drugs against SARS-CoV-2. There are no conicts to declare. COVID-19 coronavirus pandemic The PyMOL molecular graphics system The avonoids: advances in research Britannica Book of the Year Gaussian (Revision A.03) Discovery Studio Visualizer Open Source Drug Discovery and A. Lynn Frontier orbitals and organic chemical reactions This research was funded by the National Natural Science Foundation of China (grants No. 10904111, 11604238), the Tianjin Natural Science Foundation (11JCYBJC14500), and the Science and Technology Development Fund of Tianjin Education Commission for Higher Education (2019KJ175). We thank all colleagues in our laboratory for helpful discussions and technical assistance.