In Silico Prediction of Novel Inhibitors of SARS-CoV-2 Main Protease through Structure-Based Virtual Screening and Molecular Dynamic Simulation

The unprecedented pandemic of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is threatening global health. SARS-CoV-2 has caused severe disease with significant mortality since December 2019. The enzyme chymotrypsin-like protease (3CLpro) or main protease (Mpro) of the virus is considered to be a promising drug target due to its crucial role in viral replication and its genomic dissimilarity to human proteases. In this study, we implemented a structure-based virtual screening (VS) protocol in search of compounds that could inhibit the viral Mpro. A library of >eight hundred compounds was screened by molecular docking into multiple structures of Mpro, and the result was analyzed by consensus strategy. Those compounds that were ranked mutually in the ‘Top-100’ position in at least 50% of the structures were selected and their analogous binding modes predicted simultaneously in all the structures were considered as bioactive poses. Subsequently, based on the predicted physiological and pharmacokinetic behavior and interaction analysis, eleven compounds were identified as ‘Hits’ against SARS-CoV-2 Mpro. Those eleven compounds, along with the apo form of Mpro and one reference inhibitor (X77), were subjected to molecular dynamic simulation to explore the ligand-induced structural and dynamic behavior of Mpro. The MM-GBSA calculations reflect that eight out of eleven compounds specifically possess high to good binding affinities for Mpro. This study provides valuable insights to design more potent and selective inhibitors of SARS-CoV-2 Mpro.


Introduction
The current global pandemic, so called COVID-19 (COronaVIrus Disease 2019), has spread rapidly since it initially emerged in Wuhan in China, in late 2019 [1][2][3][4]. The virus called 'severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)' is responsible for the outbreak of this pandemic [5]. SARS-CoV-2 belongs to the β coronavirus subgroup of the Coronaviridae family and was found to be related to acute respiratory syndrome coronavirus (SARS-CoV) [6], which previously emerged in China in February 2003 and caused an outbreak in China and spread to several other countries [5,6]. SARS-CoV-2 specifically infects humans by causing an atypical pneumonia, which possesses specific mild to severe symptoms including dry cough, fatigue, fever, shortness of breath, severe progressive pneumonia, multiorgan failure, and eventually death [1]. The World Health Organization (WHO) has declared a state of global health emergency since the outbreak of SARS-CoV-2. According to the World Health Organization (WHO) Coronavirus (COVID-19) dashboard (https://covid19.who.int/, accessed on 30 June 2021), there have been 181,344,224 confirmed cases of COVID-19 globally, including 3,934,252 deaths worldwide, us an opportunity to explore more diverse chemical scaffolds by enhanced sampling [32][33][34][35][36]. The three-dimensional (3D-) structure of M pro is depicted in Figure 1a. Herein, we have applied target-specific virtual screening of our in-house compound database with the aim to obtain structurally diverse and potential inhibitors of SARS-CoV-2. Several compounds were identified with high inhibitory potential for M pro , and subsequently, could be tested as a treatment against COVID-19.

Validation of Docking Method by Re-docking and Cross-Docking
Prior to the virtual screening of our in-house database, re-docking and cross-docking of co-crystallized ligands were performed in order to scrutinize the efficiency of the docking method and to select the most appropriate protein file for virtual screening. The redocking results of twenty protein-ligand complexes showed that 50% of ligands were redocked with RMSD values of 0.29-1.94Å, while 30% of ligands were re-docked with RMSD ≤ 3Å. However, only four ligands showed RMSD in the range of 4 to ≥7Å. Therefore, 80% of ligands were correctly re-docked in their X-ray-determined conformations. Thus, the docking method was found efficient in predicting the experimentally determined orientations of compounds. RMSD ≤ 3.0 Å is usually considered satisfactory in redocking experiments; therefore, the results are acceptable. The re-docking results are shown in Table S1.
The cross-docking results (Table S2) showed that 40% of the ligands (ligands in 6Y2F, 6WTK, 6W79, 7BQY, 6ZRT, 7JU7, 6LU7, and 6W63) were correctly ranked between first and third position when docked in their cognate proteins, while two ligands (ligands contained in complexes with PBD codes: 7BRR and 6WNP) were ranked at fifth and seventh position in their X-ray structure. This indicates that MOE accurately ranked 50% of the ligands in good position; therefore, MOE was used in structure-based virtual screening (SB-VS) of our in-house database. The cross-docking results showed that eight proteins

Validation of Docking Method by Re-Docking and Cross-Docking
Prior to the virtual screening of our in-house database, re-docking and cross-docking of co-crystallized ligands were performed in order to scrutinize the efficiency of the docking method and to select the most appropriate protein file for virtual screening. The re-docking results of twenty protein-ligand complexes showed that 50% of ligands were re-docked with RMSD values of 0.29-1.94 Å, while 30% of ligands were re-docked with RMSD ≤ 3 Å. However, only four ligands showed RMSD in the range of 4 to ≥7 Å. Therefore, 80% of ligands were correctly re-docked in their X-ray-determined conformations. Thus, the docking method was found efficient in predicting the experimentally determined orientations of compounds. RMSD ≤ 3.0 Å is usually considered satisfactory in re-docking experiments; therefore, the results are acceptable. The re-docking results are shown in Table S1.
The cross-docking results (Table S2) showed that 40% of the ligands (ligands in 6Y2F, 6WTK, 6W79, 7BQY, 6ZRT, 7JU7, 6LU7, and 6W63) were correctly ranked between first and third position when docked in their cognate proteins, while two ligands (ligands contained in complexes with PBD codes: 7BRR and 6WNP) were ranked at fifth and seventh position in their X-ray structure. This indicates that MOE accurately ranked 50% of the ligands in good position; therefore, MOE was used in structure-based virtual screening (SB-VS) of our in-house database. The cross-docking results showed that eight proteins (PDB codes: 6Y2F, 6WTK, 6W79, 7BQY, 6ZRT, 7JU7, 6LU7, and 6W63) are appropriate for docking studies; therefore, those proteins were used in the virtual screening experiment.

Analysis of Virtual Screening Accuracy
The predictive accuracy of virtual screening was scrutinized by the ranking or the enrichment of known inhibitors (KIs, embedded in the in-house dataset) at the top-ranking position of docked libraries (Table S3). The result was examined by analyzing the percent enrichment factor (%EF) and receiver operating characteristic curves (ROC curves). These matrices are widely used to compare virtual screening results. The results showed that MOE successfully identified KIs in 6W79, 7BQY, 6ZRT, 7JU7, and 6W63 with %EF in the range of 33% to 73% at a top-100 position (Table S3), whereas 6W79 showed %EF of 53% at a top-50 position. Moreover, the ROC curve shows an AUC of 0.79-0.84 for 7JU7, 6Y2F, 6WTK, 6LU7, 7BQY, 6W63, and 6ZRT, and 0.90 for 6W79. The %EF and AUC of 6W79 were the highest among all the selected proteins. The ROC curve is displayed in Figure S1.

Pharmacokinetic Analysis
The physicochemical properties of the selected compounds showed that the molecular weight of compounds is in the range of 290 to 635 g/mol. A total of 19/30 compounds possess ≤5 rotatable bonds (RBs), while 11/30 compounds possess 6-10 RBs in their structures. The compounds have 3-13 hydrogen bond acceptor atoms (HBA) and 0-9 hydrogen bond donor atoms (HBD). The molar refractivity (MR) and topological polar surface area (TPSA) of these compounds are in the range of 74.33-169.9 and 62.32-226.83 Å 2 , respectively. These results were compared with the physicochemical properties of selected KIs. Those KIs possess 1-8 HBA atoms, 0-7 HBD atoms, and 1-22 RBs, while 4/15 KIs possess a molecular weight in the range of 549 to >681. Similarly, the TPSA of KIs was found to be in the range of 20 to >197 Å 2 . The KIs including remdesivir, lopinavir and ritonavir also have molecular weight > 600, RB = 15-23 and TPSA in the range of 120 to 203 Å 2 . According to Veber's rule of drug-likeness [39], TPSA and the number of RBs discriminate between orally active compounds and those that are not orally active for a large dataset of compounds in rats [40]. Therefore, compounds with ≤10 RBs and TPSA ≤ 140 Å 2 are predicted to have good oral bioavailability [40], while the Ghose filter further improves the predictions of drug-likeness by the following rules: the partition coefficient (LogP) of the compound should be in the range of −0.4 to +5.6, MR = 40 to 130, molecular weight = 180 to 480, and number of atoms from 20 to 70 (including HBDs and HBAs). The predicted partition coefficient (LogP octanol/water) of the selected (30)  The drug-likeness properties of selected hits were calculated based on the Lipinski rule of five [41] and Ghose [42], Veber's [39], Egan's [43], and Muegge's rules [44]. The compounds 1, 6, 12-14, 16-18, 26, 28 and 30 followed all the drug-likeness criteria given by Lipinski, while compounds 1, 4-6, 9, 12-14, 16-18, and 28 fulfilled the Ghose rules of druglikeness. Similarly, compounds 1, 6, 12-14 ). This shows that the selected hits possess comparable drug-likeness with remdesivir. Usually, substrates of biological transporters or natural products do not follow the above-mentioned rules of drug-likeness [45]. Moreover, recently several molecules were approved by the FDA in 2020. Among those approved drugs, several compounds fail on one or the other drug-likeness pharmacokinetic principle, and do not obey Lipinski, Ghose, Veber, Egan, and Muegge's filters, although this does not question the approval of these molecules. Therefore, it is critical to first look for a potent molecule, and once potency is validated, then to look for improved kinetics [46].
The bioavailability score of the compounds was in the range of 0.17-0.56, indicating moderate bioavailability. Among all the selected hits, only a few compounds (1, 2, 4, 7, 9, 10, 12, 14) showed few PAINS alerts, whereas the rest of the compounds did not show any PAINS alerts. Moreover, compounds 1, 6, 12-14, 16-18 passed the lead-likeness criteria, while compounds 2, 3-5, 7-11, 15, and 19-30 displayed few violations (i.e., MW > 300, rotors > 7, XlogP 3 > 3.5). The calculated synthetic accessibility of the compounds was in the range of 2.70 to 6.30, which reflects that these compounds are synthesizable. The bioavailability score, lead-likeness, and synthetic accessibility of compounds were compared with remdesivir, which showed that the compounds possess comparable scores. The bioavailability score of remdesivir is also 0.17 and synthetic accessibility = 6.33, and remdesivir depicted two violations in lead-likeness (i.e., MW > 350, Rotors > 7). The predicted physiological properties, pharmacokinetic profiles, drug-likeness, and medicinal properties of the selected compounds are tabulated in Tables S5-S9.
plex. Moreover, the stability of the M -12 complex was also affected due to a continuous increase in the RMSD after 40 ns. In contrast, the M pro -13 complex mediated friction between 30 and 40 ns; however, it showed minimal effect on the system's stability. Similarly, the RMSD of the M pro -17 complex was stable up to 50 ns; however, it increased at 50-60 and 75-90 ns. We observed that the M pro -8 complex depicted a rapid increase in the RMSD between 20 and 90 ns (and therefore destabilized the system); however, after a drastic increase, the RMSD was stabilized after 90 ns. The M pro -28 complex remained relaxed until 55 ns; however, the RMSD was increased after 55 ns and remained elevated throughout the simulation. The sudden increase indicates the fluctuation in the stability of the complex. Altogether, the results indicate that M pro -3, M pro -8, M pro -11, M pro -18, and M pro -28 complexes attained more variation as compared to the apo-M pro , while M pro -8, M pro -11, and M pro -28 complexes were found unstable till the end of simulation and reached a maximum RMSD of 7Å, 4Å, and 3.9Å, respectively, as compared to apo-M pro and M pro -X77 complex. During the simulation, no destruction in the simulated complexes (both apo and ligand-bound forms) was observed, which confirms the significance of the simulation. The RMSD graphs of all the complexes are shown in Figure 2. 8, M -11, M -18, and M -28 complexes revealed a similar energy pattern (in the range of -5400 kcal/mol to −5200 Kcal/mol), whereas the M pro -10 and M pro -13 complexes possess slightly higher total energy (between −5500 kcal/mol and −5700 kcal/mol) than the apo-M pro . The M pro -1, M pro -3, M pro -6, M pro -12, and M pro -17 complexes showed increased total energy ranging from −5600 kcal/mol to −5800 kcal/mol ( Figure 4). The inhibited states of M pro shared similar patterns and variations as compared to the apo form, therefore showing the effects of deviation on the structure of the protein created by each inhibitor.  The total energy of the apo-M pro was stable (energy = −5600 kcal/mol), while the M pro -X77 exhibited slightly lower energy (−5400 kcal/mol) as compared to apo-M pro . The M pro -8, M pro -11, M pro -18, and M pro -28 complexes revealed a similar energy pattern (in the range of −5400 kcal/mol to −5200 Kcal/mol), whereas the M pro -10 and M pro -13 complexes possess slightly higher total energy (between −5500 kcal/mol and −5700 kcal/mol) than the apo-M pro . The M pro -1, M pro -3, M pro -6, M pro -12, and M pro -17 complexes showed increased total energy ranging from −5600 kcal/mol to −5800 kcal/mol ( Figure 4). The inhibited states of M pro shared similar patterns and variations as compared to the apo form, therefore showing the effects of deviation on the structure of the protein created by each inhibitor.
12, M pro -17, and M pro -18 complexes depicted lower flexibility due to the differential dynamics upon ligand binding. The secondary structures with loops are responsible for the fluctuation in the RMSF at different levels, justifying the residual flexibility. The flexibility of complexes with the selected eleven hits varies as compared to the apo-M pro and X77inhibited states ( Figure 3).
The total energy of the apo-M pro was stable (energy = −5600 kcal/mol), while the M pro -X77 exhibited slightly lower energy (−5400 kcal/mol) as compared to apo-M pro . The M pro -8, M pro -11, M pro -18, and M pro -28 complexes revealed a similar energy pattern (in the range of -5400 kcal/mol to −5200 Kcal/mol), whereas the M pro -10 and M pro -13 complexes possess slightly higher total energy (between −5500 kcal/mol and −5700 kcal/mol) than the apo-M pro . The M pro -1, M pro -3, M pro -6, M pro -12, and M pro -17 complexes showed increased total energy ranging from −5600 kcal/mol to −5800 kcal/mol ( Figure 4). The inhibited states of M pro shared similar patterns and variations as compared to the apo form, therefore showing the effects of deviation on the structure of the protein created by each inhibitor.

Protein Motions and Trajectories Clustering
The dynamic impact of eleven hits on the structure of M pro is shown in Figure 5. The structural changes in each complex due to the protein-ligand binding was observed through principal component analysis (PCA). The significant dominant motions ( Figure 5) are shown in the first three eigenvectors, while the others indicated localized fluctuation. In the apo-M pro , a total of 48% of variances were contributed by the first three eigenvectors to the total observed motion. Unlike the apo-M pro , the inhibited states showed different behavior of motion. In inhibited states, compounds 8 and 12 showed 60%, compounds 18, 1 and X77 reflected 57% and 55%, respectively, whereas compounds 6, 28, and 11 showed 52-50% of total motion. The total motion of compounds 3 and 17 was 40%, while compounds 10 and 13 demonstrated least motion of 38% and 26%, respectively. These structural behavior clearly demonstrated the structural rearrangement of the protein upon ligand binding.
vectors to the total observed motion. Unlike the apo-M pro , the inhibited states showed different behavior of motion. In inhibited states, compounds 8 and 12 showed 60%, compounds 18, 1 and X77 reflected 57% and 55%, respectively, whereas compounds 6, 28, and 11 showed 52-50% of total motion. The total motion of compounds 3 and 17 was 40%, while compounds 10 and 13 demonstrated least motion of 38% and 26%, respectively. These structural behavior clearly demonstrated the structural rearrangement of the protein upon ligand binding.
The reliability of attributed motions was achieved by plotting the two initial eigenvectors of each trajectory against each other. During the simulation production run, the flipping over conformation was shown by color blue to red. The dots represented each frame from blue to red. To understand the conformational transformation of the complexes, a 2D subspace was mapped from the trajectories using PC1 and PC2. Figure 6 clearly shows that each complex acquired two conformational states on the subspace differentiated by the colors (blue and red). The unstable conformational state (shown in blue) can be easily separated in neared convergence to obtain a stable conformational state (shown in red). Subsequently, the apo-M pro showed more energetic conformation, while the inhibited states showed stable energy conformation with different periodic jumps. The M pro -X77 complex reflected very stable lower energy conformation as compared to apo-M pro , while the rest of the inhibitors followed the same pattern and acquired stability with lower energy states.  The reliability of attributed motions was achieved by plotting the two initial eigenvectors of each trajectory against each other. During the simulation production run, the flipping over conformation was shown by color blue to red. The dots represented each frame from blue to red. To understand the conformational transformation of the complexes, a 2D subspace was mapped from the trajectories using PC1 and PC2. Figure 6 clearly shows that each complex acquired two conformational states on the subspace differentiated by the colors (blue and red). The unstable conformational state (shown in blue) can be easily separated in neared convergence to obtain a stable conformational state (shown in red). Subsequently, the apo-M pro showed more energetic conformation, while the inhibited states showed stable energy conformation with different periodic jumps. The M pro -X77 complex reflected very stable lower energy conformation as compared to apo-M pro , while the rest of the inhibitors followed the same pattern and acquired stability with lower energy states. Pharmaceuticals 2021, 14, x FOR PEER REVIEW 13 of 24

Metastable to Native State Transition Pathway
The transition states of the apo-M pro and inhibited M pro complexes were studied using the free energy landscape (FEL). The FEL plot was constructed from the first two eigenvectors of the trajectory time to explore the transition mechanism from the metastable state to the native state. The lowest energy states of each complex were examined to investigate the structural changes. The apo-M pro showed a significant change in the energy states as compared to the inhibited states (represented by red, yellow, green, and blue in Figure 7). The highest energy levels and the metastable stage in the plots are shown by red and blue colors, respectively. The apo-M pro was stable as compared to the inhibited states because the red color (high energy state) is more prominent in the inhibited states (X77, 1, 3, 6, 8, 10-13, 17-18, and 28). The compounds 3, 8, 6, 11, 18, and 28 showed the highest transition states due to the interaction with the active site domain of M pro . The apo-M pro acquired only one state with no energy barriers, and similarly, compound 13 showed a pattern like apo-M pro due to the sliding of compound 13 from the active pocket because of weak interactions. Moreover, compound 10 acquired two states with a stable energy level for the maximum time (shown in yellow). The reference ligand, X77, remained mostly in the high energy state, which confirmed the effect on the stability of the protein structure due to the rearrangement of the bonds upon binding with X77. Figure 7 depicts that the apo-M pro remains in the green and yellow energy states, while ligandinhibited complexes are found in the high energy state (red) for most of the simulation time. The inhibition of M pro by the selected hits is evident by FEL, which clearly shows the structural rearrangement of the protein upon binding with small drug-like molecules. The

Metastable to Native State Transition Pathway
The transition states of the apo-M pro and inhibited M pro complexes were studied using the free energy landscape (FEL). The FEL plot was constructed from the first two eigenvectors of the trajectory time to explore the transition mechanism from the metastable state to the native state. The lowest energy states of each complex were examined to investigate the structural changes. The apo-M pro showed a significant change in the energy states as compared to the inhibited states (represented by red, yellow, green, and blue in Figure 7). The highest energy levels and the metastable stage in the plots are shown by red and blue colors, respectively. The apo-M pro was stable as compared to the inhibited states because the red color (high energy state) is more prominent in the inhibited states (X77, 1, 3, 6, 8, 10-13, 17-18, and 28). The compounds 3, 8, 6, 11, 18, and 28 showed the highest transition states due to the interaction with the active site domain of M pro . The apo-M pro acquired only one state with no energy barriers, and similarly, compound 13 showed a pattern like apo-M pro due to the sliding of compound 13 from the active pocket because of weak interactions. Moreover, compound 10 acquired two states with a stable energy level for the maximum time (shown in yellow). The reference ligand, X77, remained mostly in the high energy state, which confirmed the effect on the stability of the protein structure due to the rearrangement of the bonds upon binding with X77. Figure 7 depicts that the apo-M pro remains in the green and yellow energy states, while ligand-inhibited complexes are found in the high energy state (red) for most of the simulation time. The inhibition of M pro by the selected hits is evident by FEL, which clearly shows the structural rearrangement of the protein upon binding with small drug-like molecules. The ligand-bound complexes (inhibited states) displayed more conformational transitions as compared to the free state. Various metastable states showed conformational changes in the M pro structure in the ligand-inhibited complexes. The protein structure was ensembled at a distinct nanosecond time scale. In Figure 7, the crucial areas in the structures are shown in shaded form. The X and Y coordinates were obtained from the metastable states from all the trajectories with their respective frame number and time (ns), which are tabulated in Table 1. ligand-bound complexes (inhibited states) displayed more conformational transitions as compared to the free state. Various metastable states showed conformational changes in the M pro structure in the ligand-inhibited complexes. The protein structure was ensembled at a distinct nanosecond time scale. In Figure 7, the crucial areas in the structures are shown in shaded form. The X and Y coordinates were obtained from the metastable states from all the trajectories with their respective frame number and time (ns), which are tabulated in Table 1.

Dynamic Cross-Correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominantnegative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable states (ns) and frame number is presented in green and purple, respectively. The metastable states and the native structures of M pro are depicted in cartoon model in blue and red, respectively. The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro . 4 and energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable states (ns) and frame number is presented in green and purple, respectively. The metastable states and the native structures of M pro are depicted in cartoon model in blue and red, respectively. The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro .
5 loops showed positive correlation motion, while in the M pro -1 complex, Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. High and low energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable states (ns) and frame number is presented in green and purple, respectively. The metastable states and the native structures of M pro are depicted in cartoon model in blue and red, respectively.

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro .
2, Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. High and low energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable states (ns) and frame number is presented in green and purple, respectively. The metastable states and the native structures of M pro are depicted in cartoon model in blue and red, respectively.

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro . Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. High and low energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable states (ns) and frame number is presented in green and purple, respectively. The metastable states and the native structures of M pro are depicted in cartoon model in blue and red, respectively.

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro . Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. High and low energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable state (ns) and frame number is presented in green and purple, respectively. The metastable states and th native structures of M pro are depicted in cartoon model in blue and red, respectively.

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the func tional displacements of the protein's interactive atoms as a function of time. The apo-M pr reflected more positive correlation motion during 100 ns of simulation, while the domi nant-negative correlation motion of the loop was observed. The inhibited M pro demon strated variation in correlated motion, where maximum residues of the inhibited M pr showed positive correlation motion compared to the apo-M pro . The correlation motion o all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 dis played negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation mo tion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion while β1, β2, and β4 displayed positive correlation motion. On the other hand, com pounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5 while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. Th compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak neg ative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4 Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational change were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro .
3, Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. High and energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high en states are presented in dark blue, yellow, and red, respectively. The time of the metastable s (ns) and frame number is presented in green and purple, respectively. The metastable states an native structures of M pro are depicted in cartoon model in blue and red, respectively.

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the f tional displacements of the protein's interactive atoms as a function of time. The aporeflected more positive correlation motion during 100 ns of simulation, while the d nant-negative correlation motion of the loop was observed. The inhibited M pro dem strated variation in correlated motion, where maximum residues of the inhibited showed positive correlation motion compared to the apo-M pro . The correlation motio all the systems is graphically presented in Figure 8. The overall motions in each sys are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 played negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation tion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation mo while β1, β2, and β4 displayed positive correlation motion. On the other hand, c pounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correla motion at β1 to β2. However, compound 17 did not show significant positive correla motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak ative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion a Hence, the internal dynamics of the protein have substitution effects upon ligand bin with the protein. The results indicate that dynamic variability and conformational cha were caused by small inhibitors, therefore revealing the affinity of ligands toward M 4, and Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. energy states are shown in distinct colors in the plot. The minimal, intermediate, and th states are presented in dark blue, yellow, and red, respectively. The time of the me (ns) and frame number is presented in green and purple, respectively. The metastable native structures of M pro are depicted in cartoon model in blue and red, respectively. The dynamic cross-correlation matrix (DCCM) was constructed to elabor tional displacements of the protein's interactive atoms as a function of time. T reflected more positive correlation motion during 100 ns of simulation, whi nant-negative correlation motion of the loop was observed. The inhibited M strated variation in correlated motion, where maximum residues of the in showed positive correlation motion compared to the apo-M pro . The correlatio all the systems is graphically presented in Figure 8. The overall motions in are dominated by the correlated motions. In the X77-inhibited M pro , the β1 played negative correlation motion and ϒ4 and ϒ5 loops showed positive cor tion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correla while β1, β2, and β4 displayed positive correlation motion. On the other pounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, while the β1, β2, β3, and β4 regions acquired apparent positive correlation compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pr displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negativ motion at β1 to β2. However, compound 17 did not show significant positiv motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed ative correlation motion in β1 and ϒ1 regions, and slight positive correlation m Hence, the internal dynamics of the protein have substitution effects upon lig with the protein. The results indicate that dynamic variability and conformatio were caused by small inhibitors, therefore revealing the affinity of ligands to 5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro .

and
Pharmaceuticals 2021, 14, x FOR PEER REVIEW 15 of 24 Figure 7. The free energy landscape (FEL) of free state and inhibited states is shown. High and low energy states are shown in distinct colors in the plot. The minimal, intermediate, and the high energy states are presented in dark blue, yellow, and red, respectively. The time of the metastable states (ns) and frame number is presented in green and purple, respectively. The metastable states and the native structures of M pro are depicted in cartoon model in blue and red, respectively.

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4. Hence, the internal dynamics of the protein have substitution effects upon ligand binding with the protein. The results indicate that dynamic variability and conformational changes were caused by small inhibitors, therefore revealing the affinity of ligands toward M pro .
5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in Pharmaceuticals 2021, 14, x FOR PEER REVIEW Figure 7. The free energy landscape (FEL) of free state and in energy states are shown in distinct colors in the plot. The minim states are presented in dark blue, yellow, and red, respective (ns) and frame number is presented in green and purple, respe native structures of M pro are depicted in cartoon model in blu

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) was tional displacements of the protein's interactive atoms a reflected more positive correlation motion during 100 nant-negative correlation motion of the loop was obse strated variation in correlated motion, where maximu showed positive correlation motion compared to the ap all the systems is graphically presented in Figure 8. Th are dominated by the correlated motions. In the X77-i played negative correlation motion and ϒ4 and ϒ5 loop tion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 refle while β1, β2, and β4 displayed positive correlation m pounds 8, 11, and 18 showed negative high correlation while the β1, β2, β3, and β4 regions acquired apparen compounds 6, 10, and 13 showed similar patterns as com displayed insignificant positive correlation motion in ϒ motion at β1 to β2. However, compound 17 did not sho motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, c ative correlation motion in β1 and ϒ1 regions, and slight Hence, the internal dynamics of the protein have substit with the protein. The results indicate that dynamic varia were caused by small inhibitors, therefore revealing the

Dynamic Cross-correlated Map Analysis
The dynamic cross-correlation matrix (DCCM) tional displacements of the protein's interactive atom reflected more positive correlation motion during 1 nant-negative correlation motion of the loop was strated variation in correlated motion, where max showed positive correlation motion compared to th all the systems is graphically presented in Figure 8 are dominated by the correlated motions. In the X played negative correlation motion and ϒ4 and ϒ5 l tion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 r while β1, β2, and β4 displayed positive correlatio pounds 8, 11, and 18 showed negative high correla while the β1, β2, β3, and β4 regions acquired appa compounds 6, 10, and 13 showed similar patterns as displayed insignificant positive correlation motion motion at β1 to β2. However, compound 17 did no motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermo ative correlation motion in β1 and ϒ1 regions, and sl Hence, the internal dynamics of the protein have su with the protein. The results indicate that dynamic v were caused by small inhibitors, therefore revealing  The dynamic cross-correlation matrix (DCCM) was constructed to elaborate the functional displacements of the protein's interactive atoms as a function of time. The apo-M pro reflected more positive correlation motion during 100 ns of simulation, while the dominant-negative correlation motion of the loop was observed. The inhibited M pro demonstrated variation in correlated motion, where maximum residues of the inhibited M pro showed positive correlation motion compared to the apo-M pro . The correlation motion of all the systems is graphically presented in Figure 8. The overall motions in each system are dominated by the correlated motions. In the X77-inhibited M pro , the β1 and β2 displayed negative correlation motion and ϒ4 and ϒ5 loops showed positive correlation motion, while in the M pro -1 complex, ϒ2, ϒ3, and ϒ4 reflected negative correlation motion, while β1, β2, and β4 displayed positive correlation motion. On the other hand, compounds 8, 11, and 18 showed negative high correlation motion in loops ϒ3, ϒ4, and ϒ5, while the β1, β2, β3, and β4 regions acquired apparent positive correlation motion. The compounds 6, 10, and 13 showed similar patterns as compared to the apo-M pro , while they displayed insignificant positive correlation motion in ϒ4 and ϒ5 and negative correlation motion at β1 to β2. However, compound 17 did not show significant positive correlation motion in ϒ3, ϒ4, and β4 regions of M pro . Furthermore, compound 28 showed a weak negative correlation motion in β1 and ϒ1 regions, and slight positive correlation motion at β4.

Binding Free Energy Calculations
The binding free energy of each ligand was estimated to quantitatively compare the energy differences of the selected hits (from the in-house database) with X77. The binding free energy was computed from the last 1000 frames of the 100 ns of MD trajectory. MM-GBSA analysis was performed for each system by calculating each contributing energy, such as van der Waals (∆VDW), total electrostatic energy (∆EET), polar and non-polar contributions (∆EGB), and non-polar solvation energy (SASA) ( Table 2). The MM-GBSA results showed variation in energies among X77 and the eleven molecules. The effect is high in terms of total and electrostatic energies. The reference inhibitor, X77, exhibited

Binding Free Energy Calculations
The binding free energy of each ligand was estimated to quantitatively compare the energy differences of the selected hits (from the in-house database) with X77. The binding free energy was computed from the last 1000 frames of the 100 ns of MD trajectory. MM-GBSA analysis was performed for each system by calculating each contributing energy, such as van der Waals (∆VDW), total electrostatic energy (∆EET), polar and non-polar contributions (∆EGB), and non-polar solvation energy (SASA) ( Table 2). The MM-GBSA results showed variation in energies among X77 and the eleven molecules. The effect is high in terms of total and electrostatic energies.  . Compounds 1, 3, and 17 also showed appropriate binding potential with M pro by making stable complexes with ∆G TOTAL of −22.7848 Kcal/mol, −22.9067 Kcal/mol, and −22.0295 Kcal/mol, respectively. However, compounds 10 and 13 reflected the lowest total binding free energies (compound 10 = −6.4968 Kcal/mol and 13 = −9.4012 Kcal/mol) due to poor binding interactions within the active site of M pro .
in the library are given in SMILE format in the Supplementary Materials. Moreover, fifteen known inhibitors (KIs) were also added in our in-house database as positive controls (Table S11). The 3D-structure of each ligand (in mol2 format) was imported into MOE compound database, where Wash module of MOE was used to add hydrogen atoms and partial charges (based on Amber12:EHT force field) on each structure, and the structures were minimized with the same parameters as discussed above.

Structure-Based Screening by Molecular Docking
After the preparation of protein and ligand files, molecular docking was performed by Triangle Matcher docking algorithm and London dG scoring function [48,51,52]. The active site/ligand binding site was defined on the co-crystallized ligand in each protein. For redocking, the ligands were extracted from each protein and re-docked in its cognate binding protein (with the above-mentioned settings), and the results were quantified by calculating root mean square deviation (RMSD) between the docked and X-ray conformation of each ligand. Similarly, cross-docking was performed by docking all the twenty (extracted) ligands in each of twenty proteins and results were examined by ranking (at top position) of compounds in their native X-ray crystal structure. By default, thirty docked conformations of each ligand were obtained. The virtual screening of in-house database was performed on those PDB files that displayed good results in cross-docking experiment.

Analysis Measures and Conformational Sampling after Virtual Screening
The inhibitor with the most potential against SARS-CoV-2 M pro was chosen after virtual screening by consensus approach. The in-house library was docked in eight protein structures individually. Later, each docked library was sorted based on the docking score, and those compounds that were ranked mutually in 'Top-100' position in at least 50% of the structures were declared as potential 'Hits'. The optimal binding modes of the selected compounds were chosen through conformational sampling. The docked orientation of each compound found analogous in all the proteins was considered as the possible binding mode. The interactions of ligand were visualized by Protein-Ligand Interaction Fingerprints (PLIF) [48] of MOE, which calculates several types of interactions between protein and ligands including H-bonds, water-mediated protein-ligand bridging, ionic interactions, surface contacts, metal ligation, and arene attraction in 2D format.

Prediction of Pharmacokinetic Properties
After virtual screening, the pharmacokinetic (ADMET: absorption, distribution, metabolism, excretion, and toxicity) behavior of the selected compounds was studied through Swis-sADME [37] and admetSAR [38], which predicts ADMET properties and drug-likeness of small molecules by using physicochemical descriptors.

Molecular Dynamic Simulation
The atomic coordinates of PDB ID: 6W79 [47] were chosen for the molecular dynamic simulation studies. Thirteen systems were generated for MD simulation, including apo form of 6W79 (apo-M pro ), 6W79 in complex with co-crystallized ligand (M pro -X77), and 6W79 in complex with docked conformations of eleven hits. The apo-M pro and M pro -X77 complex were used as positive controls. The possible overlaps/clashes in the initial structure were eliminated by minimizing the structure with 10,000 cycles of steepest descent [53] (macromolecule was frozen), followed by 20,000 cycles of conjugate gradient method [54]. LEaP module of AMBER20 [55] was used to add the missing hydrogen atoms. To keep the systems neutral, counter-ions from OPC model [56] were added. A truncated octahedral box of the OPC water model [57] was added to all the systems with a 10 Å buffer (8 Å cut-off was used to compute the pairwise interactions, the van der Waals, and direct Coulomb forces). Long-range electrostatic forces were treated with the particle mesh Ewald (PME) algorithm [58]. The intermolecular interactions were calculated by ff19SB [59]. In preparation runs, Langevin thermostat [60] was used with 1 ps −1 friction constant, while Berendsen thermostat [61] was used in the production runs. MD simulation was accelerated using the PMEMD CUDA version in GPU cores. Before running MD production, all the systems were heated for 400 ps, followed by equilibration of up to 2000 ps in the NVT ensemble at 300 K. The conditions applied in the simulation of all systems are given in Table 3.

Post-Dynamic Evaluation
The coordinates of all the simulated systems were extracted from the generated trajectories after every 1 ps and analyzed by PTRAJ [62] module of the AMBER20. The Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), and radius of gyration (Rg) of all the systems were calculated by CPPTRAJ module of AMBER20 on Cα atoms via Equations (1)-(3).
In RMSD calculation, N = the number of atoms, mi = the mass of atoms, Xi = the target atom i vector coordinate, Yi = the reference atom i vector coordinate, and M = the total mass.
The RMSF of selected atom i was calculated as: the atomic positions averages over the total input frames (denoted by x).
The Rg of N number of atoms was calculated: the atomic position was denoted by ri, and the mean position was denoted by r m of all the atoms. The Altona and Sundaralingam method [63] was used to calculate the five-membered ring pucker. Standard deviation and averages were reported in the analysis utilities, with proper cyclic averages being computed for periodic values (torsions). Furthermore, the total energy of all the systems (apo-and inhibited states) was calculated.

MD Trajectories Unsupervised Clustering and Free Energy Landscape
Principal component analysis (PCA, focuses on matrix covariance) was used to demonstrate atom movement and protein loop dynamics. The internal motions of the systems were analyzed by PCA approach of CPPTRAJ. The atomic coordinates of eigenvectors and the positional covariance matrix were calculated. The orthogonal coordinate transformation was used to obtain the eigenvalue diagonal matrix by the diagonalizing of the matrix. The principal components were obtained based on eigenvalues and eigenvectors to emphasize the motion of the atoms in MD simulation trajectories. The isolated first principal components, PC1 and PC2, showing the largest variation in the data, were utilized for the free energy landscape (FEL) using Equation (4) from 100 bins of the data population. The energies were calculated in kcal/mol at 300 •K.
where K B = Boltzmann's constant, T = specified temperature, N i = bin i population, and N Max = most populated bins. The artificial barrier population size of 0.5 was applied to the bins with no population.

Dynamic Cross-Correlation Analysis (DCC)
The dynamic cross-correlation map (DCCM) method was used to obtain the Cα atom's time subordinate movements caused by the attachment of a small molecule (inhibitor) with the protein. The correlation matrix was derived by observing the Cα atoms' correlated and anti-correlated motions of each system. DCCM was calculated by Equation (5).
where C ij = time correlated data between the atoms i and j in a protein. We used 0.002 ns interval to construct the matrix of Cα from the 10,000 snapshots. The positive and negative values indicate the correlated and anti-correlated motion during the MD simulation, respectively, in the matrix plot.

MM/GBSA Free Energy Calculation
In MD simulation, free energy calculations give quantitative production of proteinligand binding energies. The binding energy (G bind ) was calculated by Equation (6).
where G R+L represents the M pro in complex with inhibitors, while G R and G L represent the apo-M pro and inhibited M pro , respectively. In the generalized born surface area (MM/GBSA) approach, each free energy term in Equation (6) was calculated using Equation (7).
where E bond represents bond angles and dihedral energy, E vdw and E elec indicate the contribution of van der Waals and electrostatic energy, respectively, while the related polar and non-polar contribution of solvation energy are reported as G GB and G SA . T and S s show the absolute temperature of the system and the solute entropy, respectively. The performance of the MMGBSA algorithm is based on the specificity of the forcefield and inhibitor's partial charges, the specificity of protein-inhibitor complex, MD simulation, inner dielectric constant, and the docking pose number based on top scoring. Here, the binding free energies of each system were calculated by MM/PB(GB)SA model of GBSA. The solvent probe of 2 Å radius was used, and the radii were used to optimize the topology files.

Data Analysis
The results were analyzed by MOE [48], UCSF Chimera [64], and Pymol [65]. The average structures were extracted from structure ensembles of the lowest energy. All the analysis graphs were plotted using Origin pro [66] and GnuPlot [67].

Conclusions
The main protease or chymotrypsin-like protease of SARS-CoV-2 is considered to be a potential anti-viral drug target. We have employed an efficient structure-based virtual screening protocol to search for novel inhibitors of SARS-CoV-2. The binding potential of several compounds was tested on multiple structures of M pro , and consensus strategy was applied to select the most promising binders. Based on the physiological and pharmacokinetic behavior and protein-ligand binding pattern, eleven compounds were identified as good inhibitors of M pro . Therefore, the structural and dynamic behavior of M pro upon binding with those eleven compounds was further explored through molecular dynamic simulation. Based on the MM-GBSA calculations, two compounds (11 and 28) were retrieved with the highest binding affinities for M pro , whereas six compounds (3, 8,  18, 6, 1, and 17) showed good binding affinities for M pro . Based on our in silico findings, we suggest that these compounds can inhibit the replication of SARS-CoV-2 by specifically inhibiting its M pro enzyme. Therefore, these compounds can act as potential anti-viral candidates against SARS-CoV-2. However, further in vitro testing is required to confirm these in silico results.