Next, to ascertain the stability of the crystal ligand’s with M^pro, we conducted a redocking process centering GC-376* and the grid that was built previously. Throughout this procedure, default parameters were strictly followed, which involved assigning the van der Waals radii of nonpolar ligand atoms to 0.8, setting a partial charge cutoff of 0.15, and employing full-flexible docking. The docking operation utilized the Glide program accessible in Maestro (Friesner et al., 2004; Halgren et al., 2004), by specifying the extra-precision (XP) mode to obtain higher accuracy (Friesner et al., 2006). Additionally, the post-docking minimization was carried out on all generated poses obtained from the previous docking calculation. Subsequently, the top ten conformations for each ligand were selected for further analysis.

Following the insights gained from the redocking procedure, we proceeded to conduct the actual docking calculations for the five different Rutan compounds. For this purpose, all the preprocessed 3D conformations of the Rutan compounds, which were generated using the LigPrep calculations (740), were chosen. These selected conformations were then subjected to the docking procedure, utilizing the previously established grid and the predefined Glide parameters. The identification of the best docking hits was based on a comprehensive evaluation, considering the following criteria: (i) graphical investigation on the binding site of M^pro, (ii) binding mode comparison with crystal ligand; i.e., comparison of the binding modes of the docked ligands with that of the crystal ligand to gauge their similarity and accuracy. (iii) glide scores: the glide scores obtained from the docking output are an integral part of the docking calculations and were utilized as a quantitative measure to rank and evaluate the binding affinities of the inhibitors. (iv) bonded and non-bonded interactions: detailed scrutiny of bonded (covalent) and non-bonded (non-covalent) interactions between the ligands and the receptor was performed to gain insights into the stability and specificity of the binding. These multi-faceted assessments allowed us to identify and select the most promising docking outputs for further analysis.

The analysis of all trajectories obtained from the molecular MD simulations was conducted using the built-in modules, such as Simulation Quality Analysis and Simulation Interaction Diagram, available in Maestro. Furthermore, we calculated the overall count of intermolecular hydrogen bonds formed between the ligands, and the binding site residues of the M^pro was computed using the Interaction Count module.

For the estimation of Binding Free Energy (BFE) values pertaining to the M^pro–Rutan compounds, we utilized the Prime thermal_mmgbsa.py (Jacobson et al., 2004) script, an integral component of the Maestro suite. The script employs the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) method to compute the BFE values for the given M^pro–Rutan complexes using the OPLS3e force field and the Variable Solvent Generalized Born (VSGB) solvation model, using default parameters. The OPLS3e force field integrates the CM1A-BCC charge model, which combines Cramer-Truhlar CM1A charges with extensive Bond Charge Correction (BCC) parameters. Additionally, default dielectric constant cutoff values (1 for the solute and 80 for the solvent) were used. Throughout the computations, the program provided estimations for various energy components, encompassing H-bond interactions, van der Waals forces, Generalized Born (GB) solvation energies, Coulombic interactions, π-π stacking interactions, lipophilic interactions, and self-contact interaction terms. These calculations were conducted separately for the M^pro, Rutans, and the M^pro–Rutan complexes. ∆ G bind = E Complex   –   E Ligand   –   E Receptor

Furthermore, the ∆G values can be subdivided into E_lipophilic, E_{electrostatics}, and E_vdW interaction components,

Where E_{electrostatics} = E_{H bond} + E_coloumb + E_{GB_solvation} and E_vdW = E_vdW + E_π-π + E_self-contact and E_lipophilic.

The MM-GBSA-based energy estimation method in Maestro does not include the calculation of conformational entropy, thus, the entropic values were not determined.

Root Mean Squared Deviation (RMSD in Å) estimation for the M^pro–Rutan compounds: The structural deviations were computed for all the Cα atoms of each M^pro-Rutan complex throughout the MD simulation (Figure 5A). The RMSD graphs illustrate that all M^pro-Rutan complexes swiftly reached an equilibrium phase following an initial relaxation period, showing deviations ranging from approximately 1.75–2.65 Å from their starting structures. Specifically, the RMSD graphs for M^pro complexed with R5, R6, and R8 compounds exhibit a consistent equilibrium phase, demonstrating a stable equilibration phase throughout the simulation, maintaining the RMSD value of ∼1.75 Å, each. This indicates that when R5, R6, and R8 compounds bind to M^pro, they do not induce substantial conformational changes compared to their initial states. Similarly, the rmsd values for the M^pro-R7 complex show a stabilized state only after experiencing an initial fluctuation phase for ∼60 ns, later maintaining an overall RMSD value of ∼2.5 Å. However, the RMSD plot for the M^pro-R7’ complex indicates a conformational deviation occurring during earlier in the simulation period, from 0 to 170 ns, followed by reaching an equilibrium phase, sustaining rmsd values between ∼1.8 and 2.3 Å. Overall, the RMSD analysis across the five distinct simulation outputs involving M^pro-Rutan complexes suggests that the binding of Rutans to M^pro does not significantly alter the overall protein conformation. Hence, the frames obtained from the final 200 ns of the MD simulations are considered reliable and suitable for further analysis.

Root Mean Square Fluctuations (RMSF in Å) estimation for the M^pro–Rutan complexes: The structural flexibility of all the Cα atoms present in each M^pro-Rutan complexes were computed during the MD simulation (Figure 5B). Overall, the RMSF graphs displayed multiple segments of structural stability and flexibility in Cα atoms of the M^pro structure in the same trend. i.e., the M^pro exhibited (i) elevated fluctuations were observed at the loop and termini regions, as anticipated, and (ii) distinct differences in backbone fluctuations were evident among the inhibitors, especially at the binding site region. Additionally, the residues involved in the binding site region exhibited only a lower level of fluctuations for the compounds R5 to R8. This indicated that Rutans exhibiting weak binding affinity displayed moderate or weak interaction with the binding site residues, leasing to higher fluctuation compared to other complexes. Conversely, inhibitors with strong binding affinity demonstrated tight and stable interactions with the binding site residues, resulting in reduced or less fluctuation of M^pro.

Radius of Gyrations (Rg in Å) estimation for the M^pro–Rutan complexes: The Rg analysis from the MD simulation will provide comprehensive information on the change in the overall internal structure or the compactness of the macromolecule due to ligand binding. During the simulation, the overall compactness is measured based on rms-distance of all Cα atoms of M^pro from its center-of-gravity. In this study, the Rg values were extracted for all M^pro-Rutan complexes from all the trajectories and plotted (Figure 5C). The Rg graph shows that from the starting initial time period until ∼100 ns, the M^pro experienced a minor fluctuation in all the systems. Subsequently, the M^pro maintained stability throughout the simulation until the end, with the stable Rg values between ∼40.5 and 41 nm for the systems R5-R7 and R8, respectively, whereas the M^pro-R7’ maintains its Rg values at ∼42 nm. Overall, the Rg plot displays only a relatively minor change in the structural compactness during the simulation, which signifies that the binding of Rutans did not majorly alter the internal compactness of M^pro.

One work has reported that the GC-376 inhibitor effectively occupied the binding pocket of M^pro (Ma et al., 2020). This binding pocket is encompassed by a network primarily composed of hydrophobic residues, prominently featuring L27, M49, L50, F140, L141, G143, C145, H163, H164, M165, L167 G170, V186, and A191. Moreover, their study revealed that inhibitor GC-376 engaged in interactions with C145 through bisulfite adducts, leading to the formation of covalent complexes characterized by a tetrahedral arrangement. Additionally, the GC-376 inhibitor was found to establish various other interactions, notably hydrogen bond interactions, with essential residues, including T26, S46, F140, C145, H163, H164, and electrostatic interactions with N142, E166, and Q189.

Concurrently, in addition to the initial structural investigation, we also explored the H bond interaction pattern of the Rutan compounds with the binding pocket residues of M^pro towards the complex formation. The initial graphical assessment of MD trajectories focusing on the interaction between Rutans and M^pro underscores a robust interaction within the binding pocket. This interaction is characterized by a firm establishment of both hydrophobic and polar interactions, with the involvement of additional residues similar to the native crystal structure contributing to the complexity of the binding. Moreover, for a deeper insight into the stability of this complex, we systematically tracked the total number of polar interactions (Figure 6A) between M^pro and Rutans over the course of the simulations.

The graph shows that during the simulation, all five complexes have revealed a significant number of polar interactions with the surrounding residues to maintain the M^pro-Rutan complexity (Figure 6A) that helped the compounds to stay intact within the binding pocket. In this aspect, the R5 and R8 compounds exhibited as high as (∼7–8) polar interactions, while the R6 and R7’ compounds exhibited relatively lesser polar contacts, i.e., ∼6 to 7, and the R7 compound exhibited much fewer (∼5–6) in comparison with other Rutans during the simulation. The 2D polar interaction graph (Figures 6B–F) reveals that E166 and Q189 residues of M^pro establish electrostatic interactions, while C145 and T190 contribute polar interactions with all Rutan compounds, whereas H41, H164, and Q192 establish polar interaction in most of the case, providing insights into the stability and frequency of these interactions.

In summary, these findings underscore the significant impact of polar interactions between the binding site residues of M^pro and the inhibitors on complex stability. This, in turn, results in elevated binding affinity values. Notably, complexes with higher binding affinity values exhibited a greater number of polar contacts, whereas those with lower binding affinity values displayed a reduction in such contacts. These insights are indispensable for the precise optimization and augmentation of the binding capabilities of a lead compound, particularly when targeting a receptor with significant activity initially detected during screening.

The essential role played by M^pro in the SARS-CoV2 viral replication process underscores the urgent need to unravel its mechanistic intricacies and expedite the discovery of high-affinity inhibitors specifically tailored for M^pro. To achieve this, we leverage advanced computational programs, both open-source and commercial, in association with the state-of-the-art high-performance computing systems. In our research, these computational tools facilitated rapid prediction of Binding Free Energies (BFEs) for a set of five distinct Rutan compounds, leveraging their three-dimensional molecular coordinates. This approach has the potential to streamline the design and evaluation of M^pro complexed with Rutan inhibitors, promising enhanced precision and efficacy in this process. Specifically, we extensively employed the thermal_mmgbsa.py—a built-in trajectory post-processing tool in Maestro—to predict BFE values for all five M^pro–Rutan complexes, averaging results from five replicates (Table 3; Supplementary Table S1). Utilizing 100 snapshots from the last 200 ns of the MD trajectories, the predicted BFE values for these complexes are as follows: −81.11, −93.75, −67.24, −92.11, and −104.63 kcal/mol, respectively.

Apart from determining the total ΔG_bind, the program also provides a detailed insights into various interaction components, encompassing Coulombic forces, covalent binding, hydrogen bonding, lipophilic interactions, π−π stacking, Solv_GB (Generalized Born solvation), solvation surface accessibility (Solv_SA), and van der Waals forces (Table 3). These individual interacting components hold immense potential for future applications, especially in the domain of novel drug discovery for M^pro inhibition. Moreover, these components offer an extensive understanding of the intricate binding process between Rutan and M^pro. Consequently, the predicted ΔG^bind is explored for its constituent energy components, enriching our understanding of the underlying binding mechanisms.

Among these individual energy components, Coulombic, Lipophilic, and Solv_SA terms exclusively contribute the most favorable energies for complex formation, while the hydrogen bond component provides relatively less favorable energy across all Rutan compounds. Notably, these energy trends align with the total number of hydrogen bond analyses depicted in Figure 6. Conversely, the Solv_GB and vdW energy terms contribute higher unfavorable energies, while the Covalent energy terms contribute relatively less unfavorably to the R5, R6, R7, and R8 compounds, and negligibly unfavorably to the R7’ system. Overall, the analysis of individual energy components across all simulation systems involving M^pro complexed with multiple Rutan compounds highlights the substantial energies provided by the Coulombic, Lipophilic, and Solv_SA energy terms for complex formation.

In the current study, the per-residue decomposition (PRD) analysis was conducted using the Prime thermal_mmgbsa.py script, with the same number of frames and an average of five replicates used to compute the BFE values for five different Rutan compounds when bound to M^pro (Table 4 PRD). This decomposition analysis provides profound insights by showcasing the specific key residues at the interface that predominantly contributed to the BFE values due to Rutans binding. Further examination, particularly focused on the binding site residues of M^pro, unraveled intricate variations in energy profiles within a diverse ensemble of hydrophobic and polar residues, shedding light on their distinctive roles in the binding process.

While the strong binding affinities of Rutan compounds to the SARS-CoV-2 M^pro active site are encouraging, it is important to evaluate their potential off-target effects. Polyphenolic compounds are known to interact with various proteases and enzymes, which could lead to unintended biological consequences. For instance, flavonoids such as quercetin have been reported to inhibit proteases beyond their intended targets, potentially affecting normal cellular processes (Di Petrillo et al., 2022). To mitigate this concern, future studies should involve broader in vitro screenings against a panel of proteases to assess the specificity of Rutan compounds. Additionally, in vivo toxicity studies will be essential to evaluate their pharmacokinetic and pharmacodynamic profiles. By understanding the potential off-target interactions, we can design more selective derivatives to minimize adverse effects.

Here, the key residues with the maximum contribution to the complex formation (hot spot residues) are extracted from all binding site residues by assigning −2 kcal/mol (Table 4; Figure 7). Based on thiscriteria, the residues with maximum contributions for Rutan compounds are C145, E166, and Q189.