Exploration of chalcones as 3-chymotrypsin-like protease (3CLpro) inhibitors of SARS-CoV-2 using computational approaches

The main protease 3CLpro is one of the potential targets against coronavirus. Inhibiting this enzyme leads to the interruption of viral replication. Chalcone and its derivatives were reported to possess the ability to bind to 3CLpro protease in the binding pocket. This study explored an in-house database of 269 chalcones as 3CLpro inhibitors using in silico screening models, including molecular docking, molecular dynamics simulation, binding free energy calculation, and ADME prediction. C264 and C235 stand out as the two most potential structures. The top hit compound C264 was with the Jamda score of −2.8329 and the MM/GBSA binding energy mean value of −28.23 ± 3.53 kcal/mol, which was lower than the reference ligand. Despite the lower mean binding energy (−22.07 ± 3.39 kcal/mol), in-depth analysis of binding interaction suggested C235 could be another potential candidate. Further, in vitro and in vivo experiments are required to confirm the inhibitory ability. Supplementary Information The online version contains supplementary material available at 10.1007/s11224-022-02000-3.


Introduction
The COVID-19 pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus was first reported in Wuhan City, Hubei Province, China, in late 2019 [1]. Until March 7th, 2022, there have been over 445 million confirmed cases of COVID-19, including nearly 6.0 million deaths, reported to the WHO [2]. Four primary treatments, including repositioning antiviral drugs, convalescent plasma, vaccines, and other adjuvant therapies are currently applied to deal with the pandemic. Vaccines are the most effective preventative therapy, with more than 10.7 billion doses have been administered. Although the numbers of antiviral agents have been developed in the early stage of the pandemic, molnupiravir was the first oral medicine approved by the FDA since December 2021 [3]. Paxlovid-a combination of nirmatrelvir and ritonavir by Pfizer, has recently been reported to reduce the risk of hospitalization or death by 89% compared to placebo in non-hospitalized high-risk adults and being assessed by the FDA for the emergency use authorization (EUA). Nirmatrelvir was specifically designed as a covalent 3CL pro inhibitor [4].
The main protease 3CLpro (also called 3-chymotrypsinlike protease or nsp5) is one of the potential medicine targets against coronavirus. The complex of 3CLpro and papain-like protease plays a vital role in the post-translation. Blocking 3CLpro could interrupt the viral replication [5]. 3CLpro is a dimer protease with two identical protomers. Each comprises three domains: domain I (residues 10 to 99), domain II (residues 100 to 182), and domain III (residues 198 to 303). Two crucial residues at the substrate-binding site (catalytic dyad) Thua-Phong Lam and Dac-Nhan Nguyen contributed equally to this work and are co-first authors.
Flavonoids are major natural structures with a wide range of bioactivities. The flavonoid family consisted of different subgroups such as chalcone, flavone, and flavonol [13]. Chalcone structures withdraw much interest from researchers due to their broad spectrum of biological activity. Different substituents on ring A and ring B of chalcone derivatives could potentially form many hydrogen bonds, inducing the binding of chalcone toward 3CLpro. This study aimed to screen for 3CLpro inhibitors from the in-house chalcone database with in silico models, including molecular docking, molecular dynamics simulation, binding free energy calculation, and ADME prediction. The study of Mathpal et al. shared a similar screening method with various in silico filters from 3000 compounds having chalcone scaffold in the PubChem library [12]. The dataset of 269 in-house chalcone derivatives used in our study is the one of readily synthesized structures. Therefore, further in vitro experiments on the top hit structures can be performed easily.

Preparation of chalcone structures
The 2D molecular structures of 269 in-house chalcone derivatives were created by ChemDraw 18.1 [14], converted into 3D-structure, molecular energy minimized by Molecular Operating Environment (MOE) 2019.01 software [15]. In the energy-minimized process, the forcefield AMBER10:EHT was used with 0.0001 RMS kcal/mol/A 2 gradient. The individual chalcone structure was saved in *.sdf format and then transformed to *.pdbqt format using OpenBabel 2.4.1 software [3].

Preparation of 3CL pro structure of SARS-CoV-2
Up to December 2021, there are 392 SARS-CoV-2 3CL pro structures in the Protein Data Bank using the 6LU7 FASTA sequence. The first crystal structure of the SARS-CoV-2 3CL pro was with the N3 covalent inhibitor. With the aim to determine the non-covalent inhibitor of the SARS-CoV-2 3CL pro , the crystal structure of SARS-CoV-2 3CL pro (PDB ID: 7L14) with the co-crystallized ligand (2-(3-(3-chloro-5-(cyclopropylmethoxy)phenyl)-2-oxo(2H-(1,3′-bipyridine))-5-yl)benzonitrile) (Compound 26) and the resolution of 1.80 Å (R-free, 0.204, and R-work, 0.176) was selected. The IC 50 value of Compound 26 (Cmpd 26) was 0.170 nM [16]. Screening models in this study were generated based on the interaction between Cmpd 26 and 3CL pro . After the protein structure was downloaded from the PDB website, water molecules were removed. Non-polar hydrogens were subsequently added and charged using AutoDock tools 1.5.7. Only chain A of the homodimer protein was used for virtual screening; chain B and their corresponding ligand were eliminated.

In silico screening approaches
The workflow of this study was illustrated in Fig. 1. At first, 269 in-house chalcones were docked into 3CL protease (PDB ID: 7L14) by AutoDock Vina 1.1.2 [17]. The outcome was then optimized and rescored by JamdaScorer [18]. The chalcones with the top 10 Jamda score were subjected to molecular dynamics simulation using GROMACS 2021.6. The stability of the protein-ligand complex was evaluated through RMSD and RMSF values. The interaction of protein and ligand was estimated by hydrogen bond frequency and MM/GBSA binding free energy calculation. The pharmacodynamic properties of potential candidates were predicted by ADME assessment (SwissADME) [19].

Molecular docking
The docking process was performed by AutoDock Vina 1.1.2 [17], in which the binding pocket was defined to cover the whole of the reference ligand Cmpd 26 and the interacted residues. When docking, the top 1 pose of each ligand was recorded. The molecular docking results were then optimized and rescored using JamdaScorer, as the JAMDA scoring function is a novel empirical scoring function built from easily optimizable individual terms. It showed excellent pose prediction performance in the CASF-2016 docking power benchmark [18]. The potential substances were selected based on the Jamda score and protein-ligand interaction model, including hydrogen bonds, ionic bonds, arene-arene, arene-H, arene-cationic, and van der Waals contacts.

Molecular dynamics simulation
The molecular dynamics simulations (MDs) of apoprotein, reference complex, and top 10 chalcone complexes were implemented by GROMACS 2021.6 [20]. The ligand topology was created by the Swissparam tool (http:// www. swiss param. ch) [21] in the CHARMM-27 force field and then merged with the protein's coordinate to generate the complex topology file in *.gro format. This study used a dodecahedron grid box with an edge of 10 Å to perform MDs, in which the system was solvated in the TIP3P water model and neutralized by adding Na + and Cl -. The energy minimization took place in 100 ps to eliminate the system's spatial conflicts or geometrical mismatches. The system was also equilibrated in "isothermal-isochoric" and "isothermal-isobaric" conditions, with two equilibrium phases, namely, NVT and NPT ensembles. Each equilibration was conducted in 100 ps, with a temperature increased to 300 K in the NVT step using the Berendsen thermostat and a stable 1 bar (0.987 atm) pressure in the NPT step using the Parrinello-Rahman pressurization unit. After reaching equilibrium, the system was simulated in 50 ns, at a temperature of 300 K, and a pressure of 1 bar. The simulation results were recorded frame-by-frame every 0.01 ns. Simultaneously, the output trajectory was then transformed into pdb format for further analysis, such as protein-ligand interactions using the PLIF tool in MOE 2019.01 and binding free energy estimation using MM/GBSA method integrated with the GMX_MMPBSA program [22].
The stability of the protein-ligand complex was evaluated through RMSD and RMSF values. RMSD values indicated the stability of the complex during simulation by analyzing the equilibration of MD trajectories (Eq. (1)). RMSF values indicated fluctuation of specific amino acids in the interactive region (Eq. (2)). A stable structure should have an RMSD value under 3 Å [23] and an RMSF value less than 2 Å [24]:

Binding free energy calculation
The binding free energy of the ligand-protein complex was calculated based on molecular mechanical surface area (MM/GBSA) using the GMX_MMPBSA program [22]. The free energy of each system was calculated based on vacuum molecular mechanical potential (E MM ), solvation energy (G solvation ), and entropy (ΔS) in Eq. (4). The binding free energy was the difference between the protein-ligand complex free energy and the total free energy of apoprotein and ligand in an aqueous solvent (Eq. (3)):

ADME prediction
The ADME properties of the top hit after the binding free energy estimation was assessed using the SwissADME tool available at http:// www. swiss adme. ch/ [19]. SwissADME is an online tool predicting the physicochemical properties, pharmacokinetics, drug-likeness, and accessibility to pharmaceutical chemistry with various reliable and quickly predictive models (Table S3). The BOILED-Egg (Brain Or Intestinal Estimated Permeation) model predicts whether a substance could passively absorb through the intestinal epithelium or blood-brain barrier. Meanwhile, the possibility of oral medicine was tested by five laws, including Lipinski et al. [25], Ghose et al. [26], Veber et al. [27], Egan et al. [28], and Muegge et al. [29]. The lead-likeness substances were filtered by Brenk et al. criteria [30], and these substances were allowed to interact with cytochrome P450. Since the SwissADME server did not provide any screening criteria for a "good" or "bad" compound, our criteria were set to assume that the best compound would not violate any of the assessed models.

Result and discussion
Database A total of 269 in-house chalcone derivatives were redrawn, minimized, and gathered in a database. The compound C1 was the basic chalcone structure. The remaining database consisted of different substituents on the two aromatic rings that were assumed to facilitate the stability of chalcones in the protease binding pocket via multiple hydrogen bonds and multiple hydrophobic interaction types. The entire database of these chalcones with ID, SMILE structure, and other molecular properties can be found in Table S1.

Molecular docking results
The grid box for the virtual screening model was set based on the co-crystallized ligand Cmpd 26 and the corresponding interacted residues. In short, the grid box center was set at (x = −21.376, y = −4.398, z = −28.478), with the size of 18 × 14 × 18 Å. The docking model was validated throughout  the redocking process, with the result described in Fig. 2. As can be seen, the RMSD value is 1.637 Å, which reveals that the re-docked pose was highly similar to the original crystal pose, and the interaction between the co-crystallize ligand and the protease was highly conserved. Thus, the docking model is reliable and suitable for virtual screening. All 269 chalcone structures were successfully docked into the active site of 3CL pro with negative AutoDock Vina docking scores (Table S2). After re-scoring and optimizing by JamdaScorer software [18], the Jamda score (JS) of the chalcone derivatives ranges from −1.38298 to −2.94106 JS, in which the co-crystallized Cmpd 26 took a value of −2.9843 JS. The top 10 chalcones with the lowest JS score are listed in Table 1, which could be divided into four chalcone structural framework groups, namely, phenothiazine ring, O-(2-(1-pyrrolidinyl)ethyl)phenol ring, N-benzylaniline ring, and tribenzyloxybenzene ring.
Chalcone with phenothiazine ring C250, C254, C264, and C267 were 4 chalcones containing phenothiazine rings that could bind strongly to 3CL pro (Fig. 3). The Jamda scores of these four chalcones were, respectively, −2.9411, −2.7456, −2.8330, and −2.7474. These substances were bound to 3CL pro via similar residues as the interaction between co-crystallized ligand Cmpd 26 and 3CL pro . All four chalcone derivatives generated van der Waals or π-sulfur interaction toward the catalytic residue Cys145. The main carbonyl and the −NH 4 + side chain of the C267 chalcone contributed attractive charge interactions toward Glu166. Meanwhile, another catalytic residue His41 connected with three out of four chalcones (C250, C254, C264) through π-cation interaction, similar to the binding between 3CL pro and the reference ligand. In addition,  the three chalcones C250, C254, and C264 also formed a hydrogen bond with Thr25 in the ammonium side chain. In brief, the C250 structure, whose ring B formed an additional hydrogen bond with Gln192 through a 3-methoxy group, was the most potent substance in this phenothiazine group.
The ketone group of the C251 chalcone formed three hydrogen bonds with Gly143, Ser144, and the catalytic Cys145, while Gln192 formed a hydrogen bond with fluor substituent. Meanwhile, the C200 ligand generated a hydrogen bond with Cys145 to the O atom of the ethoxy group and two hydrogen bonds with Leu141 and Ser144 in the N atom of the pyrrolidine ring.

Chalcones with N-benzylaniline ring
The N-benzylaniline category included three chalcones (C214, C233, and C235). These chalcones shared similar binding models in the docking complex with 3CL pro (Fig. 5), with Jamda scores of −2.7672, −2.7819, and −2.7421, respectively. The aromatic ring of the three mentioned chalcones formed π-sulfur interaction with the Met165, while N'-aniline and O'-phenol of the C233 and C235 formed hydrogen bonds with Cys145. Besides, the hydroxy substituent in the benzyl ring contributed additional hydrogen bonds in the three complexes. In addition, C214 was stabilized furthermore by the hydrogen bond interaction between the methoxy substituent of the ring B and Gln192.

Chalcones with tribenzyloxybenzene ring
C269 was the only chalcone in the top 10 substances with a tribenzyloxybenzene ring. The bulkiness of C269 might help this ligand bind strongly to the catalytic surface of 3CL pro , resulting in a good Jamda score of −2.9229 (Fig. 6). This ligand had a highly hydrophobic surface, so it was bound to the 3CL pro mainly via hydrophobic connection, including van der Waals, π-sulfur, and π-cation interactions with Met165, His41, Leu27, Pro168, and Met49 residues.

Protein stability
The stability of 3CL pro apoprotein and its complexes with top 10 chalcone derivatives during 50-ns simulation time was evaluated through the RMSD value (Fig. 7). Nine chalcone complexes almost reached equilibrium after 15 ns, which was earlier than apoprotein and the complex of co-crystalized ligand Cmpd 26. The complex with C214 chalcone failed to stabilize, with the RMSD kept escalating after 20 ns and did not reach equilibrium at the time frame of 50 ns. Therefore, the RMSD mean value and its standard deviation were calculated from the last 35-ns time frame (Table 2). Seven structures with RMSD mean values lower than that of the apoprotein and the complex of reference ligand include C250, C269, C264, C200, C251, C267, and C235. Chalcones C233 and C254 had higher RMSD mean values than the apoprotein and the Cmpd 26 complex. C214 had RMSD mean value higher than the apoprotein but lower than the Cmpd 26 complex. The RMSD plots of individual chalcone in 3CL pro complexes were provided in Fig. S1. The stability of the protein when folding was calculated through the solvent accessible surface areas (SASA) value. An increase of SASA value during the MDs means protein structural relaxation and instability. The mean value of protein backbone SASA is described in Table 2, while the SASA graph over time of each complex is illustrated in Fig. S2. In general, the presence of the chalcone ligands did not affect the protein stability since the whole SASA values were almost the same. Moreover, chalcones C250 and C264 even stabilize the protein, with a lower SASA value than the apoprotein at around 168 nm 2 . The protein stability in complexes was reinforced in the result of the radius of gyration (Rg) ( Table 2), in which the values of 10 chalcone complexes ranged from 22.06 to 22.34 Å. Compared to the Rg value of the apoprotein (22.09 ± 0.13 Å) and the complex of reference ligand (22.45 ± 0.17 Å), the figures for ten chalcone complexes were insignificantly different and even lower than that of the Cmpd 26 reference complex. This result reveals that the protein in complexes was stable during the 50 ns of MD simulations. Figure 8 shows the RMSF values of residues in the apoprotein and in the complexes. The study focused on the residues from Gly143 to Thr190 since these amino acids contributed to the binding site of 3CL pro . As can be seen, these residues fluctuated around 1-2 Å in both apoprotein and chalcone-protein complexes, meaning the interaction with ligands did not affect the stability of the protein. The 10 chalcones could bind to the 3CL pro in a stable during 50-ns simulation time, with a detailed RMSF plot shown in Fig. S4.

Evaluating the binding ability of top hits by long MD simulation
The MD process of the two most potential chalcones (C264 and C235) was extended to 100 ns (Fig. 9). The RMSD value of both ligands started to fluctuate stably from 40 ns, with an amplitude of less than 3 Å. However, it can be seen that the ligand C264 gained more stable fluctuation in the complex with 3CL pro when compared to C234. Moreover, the protein's RMSD value reveals that the binding between 3CL pro and C235 becomes unstable after 65 ns, with an increase in protein instability. These results were corresponding with RMSF plotting. In the binding sites, residues from 143 to 189 were stable with the figures under 2 Å, while some residues in the 180-190 region fluctuated more unstably (> 2 Å). Regarding the radius of gyration (Rg), the figures for two complexes show downward trends with variation, especially on the C264 complex in 70 ns. In conclusion, these two ligands C264 and C235 also show strong binding ability toward 3CL pro of SARS-CoV-2 after a long-lasting MD process. Thus, this study suggested these ligands would be promising lead compounds for an inhibitory effect on 3CL pro of SARS-CoV-2. Nevertheless, further in vitro and in vivo experiments are required to confirm the bioactivity of C264 and C235 toward 3CL pro shortly.

Binding free energy and protein-ligand interaction
After the MDs, the 500 frames of every complex during 50 ns were used to calculate binding free energy. Figure 10 shows the MM/GBSA ∆G bind values of reference and 10 chalcone complexes. Nine chalcones could bind to 3CL pro with binding free energy values ranging from −20.03 to −28.24 kcal/mol ( Table 3). As shown in Fig. 8, MM/ GBSA ∆G bind of the C214 chalcone kept escalating during the simulation. From the 42 ns to 46 ns, the ∆G bind was stable at 0 kcal/mol. This data suggested that the C214 was not stable in the binding site and completely free from the protease. This assumption was then confirmed by using VMD to visualize the simulation process. Regarding other substances, C264 would be the most potent inhibitor since this chalcone had ∆G bind values of −28.24 ± 3.53 kcal/mol in the complex with 3CL pro , which was even lower than that of the reference ligand Cmpd 26 (−27.64 ± 3.59 kcal/mol). Therefore, the interaction between C264 and 3CL pro would be further analyzed. The individual MM/GBSA binding free energy of 10 chalcones and Compound 26 with the 3CL pro can be found at Fig. S5. During 50-ns simulation, chalcone C264 bound stably to the most interactive residues of 3CL pro , including His41  Gln189 was the most multifunctional residue. It served as both hydrogen donor and hydrogen acceptor and created hydrophobic interaction with the ligand. During the simulation, there were some other hydrogen bonds between the side chain NH 4 + or the ketone group (C = O) of the ligand and the 3CL pro . While the ketone group received hydrogen bonds from Glu166, Gln189, Thr190, and Gln192 with the frequency of 8.00%, 2.68%, 3.46%, and 0.92% in order, the side chain NH 4 + group donored hydrogen to Ser46 and Gln189 with the frequency of 7.06% and 4.08%, respectively. There was another hydrogen bond between the catalytic Cys145 and the S heteroatom in the phenothiazine ring of this chalcone with a frequency of 0.34%. The interaction heat map, the MM/GBSA binding free energy of C264 with 3CL pro , and the frequency of hydrogen bond can be observed in Fig. 11.   Besides C264, C235 would probably be a 3CL pro inhibitor. This chalcone had the ∆G bind of only −22.07 ± 3.39 kcal/mol (Table 3). This outcome could be possibly affected by the high binding free energy of the initial 10-ns time frame. In the last 10 ns of simulation, this figure tended to decrease and reached −26.36 kcal/ mol, proving this chalcone stable inside the binding site. In addition, the C235 complex had the highest number of hydrogen bonds per frame with a mean value of 1.94 ± 0.55. The O atom of the phenyl ring is formed with the His163 by the remarkably high frequency of 94.56%, while the central N atom of the ligand also received a hydrogen bond from Glu166 with the frequency of 73.35%. Besides that, Glu166 also interacted with the ligand via surface and arene attraction contacts. The interaction profile between C235 and the protease was also analyzed using the PLIF tool. As a result, this complex witnessed similar interaction of His163, Glu166, Gln189, and Asp142 at 90.62%, 87.83%, 44.11%, and 20.96%, respectively. Although having a high contact frequency with the 3CL pro , the C235 structure was quite flexible inside the binding site. This figure can be observed in the ligand RMSD chart. All the mentioned data of the C235 complex during the 50-ns simulation can be found in Fig. 12.

ADME evaluation
The ADME properties of the top ten chalcones were predicted by the SwissADME server (Table 4). Based on the preliminary criteria, there was no chalcone that satisfies all models provided by the SwissADME server. However, based on the most traditional rule of five proposed by Lipinski et al. [25], nine chalcones could become an oral drug with respect to bioavailability. The absorption ability of the ten hits was also evaluated using the BOILED-Egg model (Fig. S6) [31]. The two top hits C264 and C251 were predicted to have the ability to penetrate the blood-brain barrier to get inside the brain. The result for the high Brenk alert point observed in Table 4 because Michael acceptor in chalcone structures was 1 of 105 alert fragments proposed by Brenk et al. [30]. After the ADME evaluation, C264 was still the most potent substance as it was predicted to inhibit only CYP2C9, and chemical optimization can be conducted to resolve this problem.

Conclusion
This study used 269 in-house chalcone derivatives to screen for potential 3CLpro inhibitors by the in silico approach. Firstly, the study used the docking method to evaluate the binding ability of these chalcones toward the main protease of SARS-CoV-2, with the docking results optimized by JamdaScorer [18]. The top 10 structures were then subjected to molecular dynamics simulation, with the stability of the protein-ligand complex evaluated through RMSD and RMSF values. After that, the binding free energy of the top hits was calculated, and the pharmacokinetics of potential substances were also predicted. As a result, the structures C264 and C235 have been identified as promised candidates to inhibit the main protease of SARS-CoV-2. These selected compounds showed strong biding ability toward 3CLpro via hydrophobic and hydrogen bond interactions, with Jamda score of −2.9229 and −2.7421, respectively. The binding free energy values of C264 and C235 were respective at −28.23 ± 3.53 kcal/ mol and −22.07 ± 3.39. These lead compounds meet the ADME requirements, but the structures must be modified in the next step. Further in vitro and in vivo experiments are also required to confirm the inhibitory ability of the two chalcone structures toward 3CLpro of SARS-CoV-2.

Availability of data and material
The data used to support the findings of this study are included within the article and also available on request from the corresponding authors. Table 4 ADME properties of the top 10 chalcones retrieved from the SwissADME server