**2. Results**

#### *2.1. Comparative Analysis of PF-07321332, α-ketoamide, Lopinavir, and Ritonavir Binding to 3CLpro*

The crystal structure of 3CLpro in complex with PF-07321332, lopinavir, and ritonavir is not ye<sup>t</sup> available. Therefore, the docked complex of PF-07321332, lopinavir, and ritonavir with 3CLpro was generated using the Molecular operating environment (MOE) software. Among the top 10 docking solutions, the best ligand conformation was selected based on the docking score and binding pose. In addition, the structure of α-ketoamide co-crystalized with 3CLpro (PDB ID: 6Y2F) was retrieved from the RCSB PDB database.

To elucidate the mechanism of binding of PF-07321332 (PF), α-ketoamide (keto), lopinavir (lop), and ritonavir (rit) with 3CLpro, all complexes were subjected to 100 ns simulations in Gromacs. To assess the stability of all four systems, root mean square deviation (RMSD) values were calculated for each ligand and protein individually and in complexed form, as presented in Figure 2. RMSD of apo-3CLpro at 45 ns indicated a deviation up to 3 Å, after which the system was stabilized and oscillated in the mean range of 2.5 Å (Figure 2A). The protein in 3CLpro–PF (Figure 2B) remained stable throughout simulation with an average RMSD of 2.7 Å, while in 3CLpro–keto, the protein fluctuated in a 30–50 ns interval. The protein in 3CLpro–rit showed RMSD rising above 3 Å after 60 ns with higher fluctuation at 90 ns. The RMSD protein graph of 3CLpro–lop and 3CLpro– rit deviated from the apo-form. The 3CLpro–rit showed sinusoidal behavior until 60 ns and higher peak in the 60–80 ns time interval. Nevertheless, the RMSD of the protein in 3CLpro–lop increased gradually from 60 ns onwards (Figure 2B). It also manifested fewer fluctuations initially than apo-3CLpro did until 40 ns, and then both systems yielded somewhat similar patterns up to 75 ns, and suddenly, the protein in 3CLpro–rit showed a remarkable deviation between 75 ns and 100 ns (Figure 2B). While assessing the RMSD of all four complexes (Figure 2C), we noticed that the 3CLpro–PF complex remained stable during the simulation. The RMSD of 3CLpro–keto complex fluctuated between 20–50 ns and then gradually stabilized. 3CLpro–rit and 3CLpro–lop complexes showed similar trend as their proteins. The RMSD graph of the ligands (Figure 2D) indicated that lopinavir and ritonavir featured greater deviations in comparison with other two ligands, whereas PF manifested the least deviations. The ligand in the 3CLpro–keto complex showed high jumps around 70 ns but quickly gained stability with a RMSD value of 1.5 Å. RMSD of the ligand in 3CLpro–lop showed increasing trend from 45 ns onwards. The ligand in 3CLpro–rit manifested a sudden increase in RMSD after 20 ns and fluctuated between 3 and 4.5 Å but attained stability later. We concluded that inhibitors PF and α-ketoamide imparted stability to the protein in their complexes.

**Figure 2.** (**A**) The root mean square deviation (RMSD) graph of the apo form of 3CLpro during the 100 ns simulation. The graph suggests that protein was stable during the whole simulation, with an average RMSD of 2.2 Å. (**B**) The RMSD graphs of all four proteins in complexes in comparison with the apo form. The protein in 3CLpro–PF and 3CLpro–keto remained stable whereas in 3CLpro–lop and 3CLpro–rit, it showed more fluctuations than in the other two 3CLpro systems. (**C**) The RMSD graph of 3CLpro complexes. The 3CLpro–rit and 3CLpro–lop complexes showed similar trend as the proteins. (**D**) RMSD graphs of the ligands in complexes with 3CLpro. The ligands PF-07321332 and α-ketoamide seemed stable, whereas both lopinavir and ritonavir in 3CLpro complexes featured noticeable fluctuations. PF—PF-07321332; keto—<sup>α</sup>ketoamide; rit—ritonavir; lop—lopinavir.

To analyze the dynamic behavior of protein residues, root mean square fluctuation (RMSF) values were assessed. As compared with the apo form, RMSF (Figure 3A) of the protein in 3CLpro–rit indicated that amino acid residues 45–75 in the spanning region fluctuated by more than 3 Å. By contrast, in 3CLpro–lop, the residues of 3CLpro between positions 45 and 75 fluctuated in the range similar to that of its apo form, i.e., 2.5 Å (Figure 3A). The rest of the protein residues of both complexes oscillated equally relative to the apoprotein. Overall, the residues of the 3CLpro in complex with the ligands PF and α-ketoamide did not undergo any fluctuation, while in 3CLpro–rit complex, the protein featured the greatest fluctuations in comparison with the apo form.

To assess the compactness of protein structure, we calculated the radius of gyration (Rg). This parameter of the protein in 3CLpro–PF (Figure 3B) remained in the steady state and close to that of the apoprotein during the simulations, which means that the folding of the protein in 3CLpro–PF complex remained undisturbed. Figure 3B illustrates that the apoprotein followed a similar trend and oscillated at ~22 Å. The Rg of the protein in 3CLpro–keto fluctuated between 20–50 ns; thereafter, it retained steady state for the rest of the simulation. Rg of 3CLpro–rit remained steady up to 50 ns, then it gradually increased. 3CLpro–lop fluctuated with the highest peak of 23 Å near time point 39 ns, indicating that protein folding was disturbed. Lopinavir and ritonavir might have led to the unfolding in the 3CLpro and thus exhibited increasing trend in Rg of the protein.

Next, the number of hydrogen bonds was calculated to evaluate the protein–ligand interactions during the simulations. The PF formed three hydrogen bonds (Figure 3C), and this remained consistent throughout the time trajectory. The α-ketoamide formed two bonds initially but their number dropped to one after 50 ns. The lopinavir initially formed three hydrogen bonds with 3CLpro, but their number gradually decreased, and the hydrogen bonds disappeared around 48 ns. Lopinavir formed two bonds in the beginning

and ritonavir formed one, but these interactions were inconsistent throughout the MD simulations (Figure 3C). Therefore, the hydrogen bonds formed by lopinavir and ritonavir with 3CLpro are inconsistent and weak as compared with the interactions of 3CLpro with PF-07321332 and α-ketoamide.

**Figure 3.** (**A**) Root mean square fluctuation graphs of apo-3CLpro and of 3CLpro complexes. In comparison with the apo-form, regions spanning amino acid residues 45–75 and 150–200 in both docked complexes showed deviations. (**B**) Radius of gyration of the 3CLpro in complexes was higher than that of its apo-form. (**C**) Number of hydrogen bonds in the four complexes. PF—PF-07321332; keto—α-ketoamide; rit—ritonavir; lop—lopinavir.

#### *2.2. Interaction Analysis of PF-07321332, α-ketoamide, Lopinavir, and Ritonavir with 3CLpro*

To evaluate the ligand–protein interactions, we computed the minimum distance between interacting atoms of the ligand and protein residues. Our analysis revealed the dynamics of the interaction of PF-07321332, α-ketoamide, lopinavir, and ritonavir with 3CLpro (Figure 4). In 3CLpro–PF, three hydrogen bonds were formed (Table 1). The distance of Cys145–O2, Glu166–O3, and Gln189–H20 bonds remained within 4 Å (Figure 4A). The catalytic dyad residue Cys145 made strong bonds with atoms O2 and H9 as they remained intact during the simulation. This showed that the PF-07321332 has high affinity towards 3CLpro. The 4-methoxylindole group interacted through a hydrogen bond with the backbone of Glu166, a key residue located in the middle of the binding site. The study of Glu166 mutation, carried out in SARS-CoV 3CL protease (96% identity to SARS-CoV-2 3CLpro), has corroborated the role of this key residue in the substrate-induced dimerization [23]. The trifluoroacetamide moiety of the PF-07321332 compound involved in the Gln189–H20 bond accommodated the cyclopropyl-proline moiety to occupy the central portion of the binding site. Subsequently, the pyrrolidone moiety rearrangemen<sup>t</sup> led to bond formation of reactive nitrile group with Cys145. According to the analysis of the entire trajectory of 3CLpro–PF complex, PF-07321332 contacted Leu141 and Gln142. Leu141, Gln142, and Cys145 were part of the oxyanion loop (residues 138–145), which is lining the glutamine binding pocket and is presumably involved in the stabilization of the tetrahedral acyl transition state [24]. The ligand in 3CLpro–keto made three hydrogen bond interactions (Glu166–O48, Cys145– O41, and Met165–O48) and one arene interaction, His41–H32 (Figure 4B). The cyclopropyl moiety facilitated ligand binding in the shallow substrate-binding pocket at the surface of the 3CLpro. The α-keto group of the inhibitor interacted with Cys145. The study reported that thiohemiketal formation occurred in a reversible reaction by the nucleophilic attack of the catalytic Cys145 on the ligand [14]. In our study, α-ketoamide occupied the substrate

binding space and its carbonyl oxygen formed hydrogen bond with the main-chain oxygen of Glu166. All these interactions kept fluctuating as the distance between the interacting entities continuously changed due to the bulky nature of the ligand. Similarly, in the 3CLpro–rit complex, three hydrogen bonds and an arene-type interaction were formed (Figure 4D). The length of His41–H6 and Gly143–O3 bonds remained within 7 Å during the initial 50 ns and continuously increased onwards. In the arene-type interaction, the distance between Pro168 and the thiazolyl group fluctuated within 5–10 Å during the simulation. The hydrogen bond Thr190–S2 appeared during the initial 10 ns, and the bond length was within 5 Å, but it diminished after 10 ns (Figure 4D). Two hydrogen bonds formed in 3CLpro–lop. Initially, the bond length of His41–O38 and Glu166–O20 (Table 1) was less than 5 Å, but after 50 ns, an increase in the bond length was detected because of ligand fluctuations, as depicted in the RMSD graph (Figures 2D and 4C).

**Figure 4.** The minimum distance graph of protein–ligand interactions. The interaction distance between the protein and ligand in 3CLpro–PF (**A**) and 3CLpro–keto (**B**) remained constant. Interaction in 3CLpro–lop (**C**) initially remained constant but showed major fluctuations afterwards. In 3CLpro–rit (**D**), the lengths of all four bonds continuously increased. PF—PF-07321332; keto—α-ketoamide; rit—ritonavir; lop—lopinavir.


**Table 1.** Residue and ligand atom interactions with the bond type and energy in 3CLpro complexes.

#### *2.3. Binding Free Energy Calculation*

Protein–ligand binding affinity is essential in understanding molecular recognition. The bespoke antiviral clinical candidate PF-07321332 (NCT04756531) is a potent 3CLpro inhibitor with antiviral activity against SARS-CoV-2. The relative binding free energy of PF-07321332, α-ketoamide, lopinavir, and ritonavir with 3CLpro was computed. To calculate the binding affinity of these ligands with 3CLpro, we extracted 1000 frames from the simulation trajectories. The g\_mmpbsa (Table 2) computation revealed that the binding affinity of PF-07321332 for 3CLpro was higher in comparison with α-ketoamide, lopinavir, and ritonavir. Relative to α-ketoamide, PF-07321332 showed more considerable binding energy toward 3CLpro, almost by 28 kJ/mol, and three times more than lopinavir and ritonavir. These data showed that PF-07321332 had a stronger binding affinity for 3CLpro.

**Table 2.** Binding energy of 3CLpro and HIV-I protease complexes as computed via the g\_mmpbsa method.


## **3. Discussion**

Combating the COVID-19 pandemic requires both prevention, including vaccinebased prevention, and targeted treatment for those who ge<sup>t</sup> infected. Given the way SARS-CoV-2 is evolving into new genetic variants and the worldwide impact of COVID-19, it is likely that access to treatment alternatives is crucial, now and beyond the pandemic. The frequent mutations aid the virus to infect different hosts and increase the virulence and transmissibility of SARS-CoV-2. The RNA viruses, including SARS-CoV-2, develop one mutation per replication cycle, which means RNA viruses may adapt to new environment, but these viruses are also restricted in their ability to change their genomes because they must keep their mutation rate low in order to survive [25]. Several mutations have been reported in the viral genome so far; 106 major amino acid substitutions were found to be responsible for increased virulence and transmissibility of the virus [26]. Most of these missense mutations occurred in structural proteins of SARS-CoV-2, such as spike (S) [27] and nucleocapsid (N) proteins. Only one point mutation, G15S, has been predicted using in silico studies in 3CLpro, but this residue is far from the active site and has a destabilizing effect on the protein structure [26]. Therefore, a virus comprising such mutations with a destabilizing effect would not last for long and would have reduced virulence. Only stabilizing mutations, such as in the S protein, increase the virulence and transmissibility of the virus [28]. A recent study has reported that the significant key mutations are located at the inverted repeat regions and CpG islands. These findings showed the role of mutations in the instability of the viral genome [29].

In this study, we have explored the binding of PF-07321332, α-ketoamide, lopinavir, and ritonavir to SARS-CoV-2 3CLpro through MD simulation and MMPBSA calculation. PF-07321332 (NCT04756531) and α-ketoamide molecules that specifically bind to and inhibit SARS-CoV-2 3CLpro could be promising alternatives to fight the pandemic [14]. The comparative binding mode analysis of PF-07321332, α-ketoamide, lopinavir, and ritonavir might provide a glimpse into designing rational drugs through the modification of the inhibitors based on the residues in the active site of the enzyme.

To determine the mechanism of binding of PF-07321332, α-ketoamide, lopinavir, and ritonavir to 3CLpro, molecular docking of these molecules was done. The best docking solutions for PF-07321332, lopinavir, and ritonavir were selected based on the docking score and the best binding pose of the ligand. To understand the binding mechanism of PF-07321332 and α-ketoamide along with the two antiretroviral drugs, lopinavir and ritonavir, to 3CLpro, MD simulations were performed for 100 ns. Backbone RMSDs of the 3CLpro–PF and 3CLpro–keto complexes were stable, whereas the 3CLpro–lop complex featured a deviation at 60–80 ns and the 3CLpro–rit complex a deviation between 75 ns and 100 ns (Figure 2B). RMSD deviation of the 3CLpro protein was also analyzed in the ligand-bound form. The RMSD of 3CLpro (Figure 2C) was found to be stabler in association with either PF-07321332 or α-ketoamide than with lopinavir or ritonavir. In contrast to the other systems, 3CLpro in complex with ritonavir featured marked deviation between 88 and 95 ns of the simulation. Both ligands PF-07321332 and α-ketoamide are stabler with 3CLpro in comparison with lopinavir and ritonavir because later two underwent marked deviation during the simulation (Figure 2D). The recent clinical trials have also showed these antiretrovirals are not efficacious as they did not significantly accelerate clinical improvements in serious COVID-19 patients [30]. Our binding energy data also showed low affinity of these drugs towards 3CLpro. To assess the impact of ligand binding on protein residues, RMSF of 3CLpro were analyzed in all four complexes. With either lopinavir or ritonavir, amino acid residues 45–65 and 145–200 showed greater fluctuations in comparison with the apo-form (Figure 3A). These fluctuating regions mainly consist of binding-site residues that are interconnected by a hydrogen bond network and are involved in catalytic dyad formation between His41 and Cys145 [31–33]. However, the RMSF of the protein in 3CLpro–PF showed similar trend as its apo-form. Compactness, which determines stability of a protein, was evaluated for 3CLpro using Rg. The compactness of 3CLpro in complex with PF-07321332 and α-ketoamide remained in steady state the same as the one of its apo-form (Figure 3B). Rg of 3CLpro in its complex with ritonavir or lopinavir was higher than that of apo-3CLpro. Nonetheless, 3CLpro compactness after PF-07321332 and α-ketoamide binding remained constant throughout the simulation.

One of the most interesting examples of MD application is the drug discovery area, where this method drives experiments [34,35]. To analyze the interactions of the four molecules with 3CLpro, the minimum distance between the atoms of interacting residues of each protein and ligand (PF-07321332, α-ketoamide, lopinavir, and ritonavir) were calculated as a function of simulation time to analyze the dynamic behavior of the interacting groups. In the interaction analysis, PF-07321332 and α-ketoamide showed stronger interactions with 3CLpro, and these interactions remained intact during the simulation because the minimum distance between the interacting groups stayed almost unchanged. The interactions, i.e., Cys145–O2 and Cys145–H9 in 3CLpro–PF (Figure 4A), remained intact during the simulation. This showed that PF-07321332 was strongly bonded to 3CLpro and capable to disrupt the His41–Cys145 catalytic dyad, which, together with the N-terminus residues 1 to 7, is thought to have a vital role in proteolytic activity [36]. The Glu166 anchor hooks the ligand tightly to the central region of the binding site, which facilitate the formation of further interactions with other residues. Gly143 of SARS-CoV-2 3CLpro was reported to be the most favorable residue to form hydrogen bonds with ligands, followed by Glu166, Cys145, and His163 [37]. The amides of Gly143, Cys145, and Ser144 form the cysteine protease's canonical "oxyanion hole". The interaction of PF-07321332 and α-ketoamide with the catalytic site residues may cause the distortion of the oxyanion hole in the reaction mechanism, and it may lead to the inhibition of 3CLpro in SARS-CoV-2. Compared with PF-07321332 and α-ketoamide, both antiretroviral drugs lopinavir and ritonavir manifested weaker interactions with 3CLpro, judging by binding energy values (Table 2). PF-07321335 interacts with Cys145, thereby disrupting His41–Cys145 catalytic dyad. The α-ketoamide has the same mechanism as reported earlier carries; a nucleophilic attack of the catalytic Cys145 on the α-keto group of the inhibitor, disrupting the His41–Cys145 catalytic dyad [14]. Our binding energy data (Table 2) revealed that PF-07321335 has higher affinity towards 3CLpro. This is in accordance with PF-07321335 showing a potent inhibition of 3CLpro (NCT04756531). Lopinavir and ritonavir showed weaker interactions because only two residues, His41 and Glu166, were found to be involved in hydrogen bonding. The minimum distance of 5 Å between these interacting entities persisted up to 45 ns; after that, it fluctuated. In contrast with PF-07321335 and α-ketoamide, lopinavir and ritonavir engaged in interactions with 3CLpro, but these interactions varied substantially. Some chemical substitutions in lopinavir and ritonavir (such as the addition of an α-keto group) can act as electrophiles that may prevent the His41–Cys145 dyad formation by a nucleophilic attack of Cys145. Although in our simulations, lopinavir and ritonavir engaged in interactions with His41 and Cys145 along with other catalytic site residues, these contacts fluctuated considerably, as evident from our interaction analysis. Coronavirus protease 3CLpro lacks a C2-symmetric pocket, through which lopinavir and ritonavir bind to HIV-I protease. The absence of a C2-symmetric pocket in coronavirus 3CLpro could be the reason why these drugs do not bind to 3CLpro, which makes them ineffective in COVID-19 patients. The binding energy data (Table 2) showed that lopinavir and ritonavir have weak affinity to 3CLpro.

Conclusively, 3CLpro is an important drug target, which cleaves the viral polyprotein at 11 different cleavage sites and is required for viral maturation. The druggability of 3CLpro has already been demonstrated in various studies. The mutations in the dominant variants of SARS-CoV-2, such as B.1.1.7 and B.1.617, with higher transmissibility rate were studied and it has been demonstrated that most of the mutations have been found in the receptor binding domains or structural proteins, but no mutation has been reported in 3CLpro so far [38]. Therefore, 3CLpro is an attractive target for anti-COVID-19 drug discovery. In recent studies, many new compounds have been reported as potential inhibitors of 3CLpro [14,39–42]. We included four of the potential inhibitors in our study, compared their binding energy and discussed their interactions with the protein to identify the potent inhibitor. We found that the compound PF-07321332, currently in clinical trials, was the best of these four compounds. In addition to these molecules, sterenin M was also reported in a recent computational study [39] to bind to the same active binding site as PF-07321332 and interact with the residues of catalytic dyad, but its total binding energy profile showed that the binding of PF-07321332 with protein was much stabler than the one of sterenin M [39]. In similar studies, the identified compounds bind to the same binding site, but the stability of complex is lower than the one of PF-07321332 [40,41]. Hence, considering the current pandemic situation, it is crucial to find some potent drug candidate with good binding affinity. The MD simulation and MMPBSA are widely used for assessing protein– ligand binding affinity, but there are a few limitations of the MD simulations, such as many biochemical processes including receptor conformational shifts relevant to drug binding occur on time scales that are much longer than those amenable to simulations. In MD simulations, each atom is assigned a fixed partial charge before simulation. However, the electron clouds in surrounding atoms shift continuously according to their environments. A dynamic and responsive representation of atomic partial charges would be more accurate.
