Investigating novel thiazolyl-indazole derivatives as scaffolds for SARS-CoV-2 MPro inhibitors

COVID-19 is a global pandemic caused by infection with the SARS-CoV-2 virus. Remdesivir, a SARS-CoV-2 RNA polymerase inhibitor, is the only drug to have received widespread approval for treatment of COVID-19. The SARS-CoV-2 main protease enzyme (MPro), essential for viral replication and transcription, remains an active target in the search for new treatments. In this study, the ability of novel thiazolyl-indazole derivatives to inhibit MPro is evaluated. These compounds were synthesized via the heterocyclization of phenacyl bromide with (R)-carvone, (R)-pulegone and (R)-menthone thiosemicarbazones. The binding affinity and binding interactions of each compound were evaluated through Schrödinger Glide docking, AMBER molecular dynamics simulations, and MM-GBSA free energy estimation, and these results were compared with similar calculations of MPro binding various 5-mer substrates (VKLQA, VKLQS, VKLQG) and a previously identified MPro tight-binder X77. From these simulations, we can see that binding is driven by residue specific interactions such as π-stacking with His41, and S/π interactions with Met49 and Met165. The compounds were also experimentally evaluated in a MPro biochemical assay and the most potent compound containing a phenylthiazole moiety inhibited protease activity with an IC50 of 92.9 μM. This suggests that the phenylthiazole scaffold is a promising candidate for the development of future MPro inhibitors.

Owing to their significance in the viral life cycle and lack of related human homologues, the SARS-CoV-2 RNA-dependent RNA polymerase (RdRp) and 3C-like main protease (3CL Pro or M Pro ) have been identified as potentially promising drug targets [11,12]. The RdRp catalyzes the synthesis of viral RNA and is the target of the nucleotide analog prodrug remdesivir [11,13]. With the Food and Drug Administration's Emergency Use Authorization of remdesivir, the RdRp is currently the only SARS-CoV-2 drug target with an approved medicinal therapy [13]. M Pro , a homodimeric enzyme characterized by a catalytic His-Cys dyad, functions to cleave polyproteins pp1a and pp1ab at the recognition sequence  [12,14,15]. This function makes M Pro a potential target of therapeutic drugs [14,16]. While the Pfizer PF-07321332 covalent M Pro inhibitor has recently shown remarkably effective clinic results, there are currently no FDA-approved non-covalent inhibitors.
Thiazole and indazole derivatives have displayed a wide array of biological function and there is significant interest in their pharmacologic applications. Many bioactive natural products, such as Vitamin B1, bacitracin, and penicillin contain thiazole ring structures [31]. Additionally, synthetic thiazole-based compounds have been shown to function as antineoplastic agents, anti-HIV drugs, antifungal agents, antiparasitic agents, anti-inflammatory agents, and antiulcer agents [41,42]. Most recently, thiazole-based inhibitors have been reported for SARS-CoV-2 [43,44]. Synthetic indazole-based compounds have exhibited function as anti-inflammatory, antiarrhythmic, antitumor, antifungal, antibacterial, and anti-HIV drugs [45]. Indazole-based drugs are also being considered as SARS-CoV-2 M Pro inhibitors [46].
Thiazole-based compounds have shown promise as SARS-CoV-2 M Pro inhibitors with best-in-class compounds demonstrating K I values of 2.2 μM and IC 50 as low as 3 μM [43,44]. A machine learning approach to identifying investigational or off-market drug targets for SARS-CoV-2 identified an indazole containing compound as one of the most promising with an estimated binding affinity of À9 kcal/mol [46]. Given the wide-ranging and beneficial pharmacological impacts of thiazole-and indazole-containing compounds, and the community's interest in these drugs as potential SARS-CoV-2 M Pro inhibitors, we seek to elucidate the usefulness of thiazolyl-indazole derivatives as potential scaffolds for further pharmaceutical development and exploration [41][42][43][44][45][46]48].
In this work, we explore the potential for compounds I-V to function as reversible, noncovalent inhibitors of M Pro through chemical property prediction, biochemical assays, docking analysis and molecular dynamics simulations. For benchmarking purposes, we also include a detailed binding comparison with the viral substrate as well as a recently reported tight-binding non-covalent inhibitor X77 (PDB code 6w63: estimated K d ¼ 0.057 μM) [26]. We have identified M Pro amino acids that consistently interact with these compounds and our simulations reveal the importance of π interactions in the mechanism of binding. We demonstrate experimentally that compound I is capable of inhibiting the activity of M Pro ; and this, taken together with our binding analyses, suggests that the phenylthiazole scaffold is a good candidate for future M Pro inhibitor drug development.

Synthesis
Compounds I and II have been previously reported and their spectral details are included in Supplemental Information (Figures S1-S6) [47]. The synthesis of the thiazolyl-indazole heterocycles IV and V is being reported for the first time in this report. Each compound was prepared individually in a one pot reaction of the corresponding natural monoterpenic ketone, thiosemicarbazide, and the phenacyl bromide in refluxing ethanol conditions (Scheme 2).
The general synthetic procedure was as follows. To an ethanolic solution (50 mL of absolute ethanol), 1 equivalent each of phenacyl bromide 1 and thiosemicarbazide 2 was added to 1 equivalent of (R)-Pulegone to produce IV, or to (R)-Menthone to produce V. The reaction mixture was heated under reflux for 1 h. After evaporating the solvent, the residue was diluted with water (10 mL) and extracted with dichloromethane. The organic layer was separated, dried on anhydrous MgSO 4 , evaporated to dryness, and then purified by silica gel column Fig. 1. Subsites of the SARS-CoV-2 M Pro with the N3 inhibitor covalently bound (PDB Code 6lu7). Subsites S1 0 , S1, S2 and S4 are displayed in green, red, cyan, and yellow, respectively. Catalytic residues His41 and Cys145 are shown in black below the surface of the subsites. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) chromatography, using hexane/ethyl acetate as eluent, to obtain the corresponding thiazolyl-indazole compounds in good yield (IV: 82%; V: 70%). The structures of the two synthesized products were confirmed based on their spectral data (Figures S7 -S12). The 1 H NMR spectra showed the appearance of the thiazolic proton H 5' as a singlet at 6.82 ppm and at 6.86 ppm for IV and V, respectively. In addition, each compound displayed peaks in the range of δ H ¼ 7.29-7.92 for IV and δ H ¼ 7.29-7.82 for V corresponding to a phenyl group resonance. The 13 C NMR spectra reveal the thiazolic carbon C 5' at 102.40 ppm and 102.79 ppm for IV and V, respectively. Further NMR spectral details can be found in the Supplemental Information.

Biochemical assay
All compounds were tested using the 3CL Protease, Untagged (SARS-CoV-2) Assay Kit (BPS Bioscience, San Diego, CA). The manufacturer's protocol was followed with slight modifications. Briefly, 5X serial dilutions of test compounds or covalent inhibitor control GC376 were prepared in water. M Pro was diluted in the Assay Buffer containing 1 mM DTT at 0.5 ng/μL. The diluted enzyme solution was pipetted in duplicate into a 96-well half area opaque plate (30 μL/well) while buffer alone was pipetted into "blank wells." Dilutions of compounds I, II, and V, GC376, or 5% DMSO in water was added to wells containing the enzyme to 1X concentrations indicated in Fig. 2, S13, and S14 and incubated at room temperature for 30 min with slow shaking. Protease Substrate was added to all wells and incubated with slow shaking for an additional 4 h. Fluorescence measurements were taken using a SpectraMax i3 plate reader (Molecular Devices) with fixed excitation (360 nm) and emission (460 nm). The endpoint values were recorded using the Softmax Pro software and all raw data was normalized with the "blank" solutions.

Structure retrieval and preparation
The structure of the SARS-CoV-2 M Pro was obtained from the Protein Data Bank (PDB code 6lu7) [51]. All waters and N3 were removed. A dimerized model of M Pro was created by aligning 6lu7 polypeptides to chains A and B of Human Coronavirus NL63 M Pro (PDB code 5gwy). Previous studies have shown that the monomer is functionally inactive due to the role N-terminal residues of the second chain play in forming the active site [52]. Schr€ odinger's Protein Preparation Wizard was used together with Prime, Epik, and PROPKA to prepare the protein [53][54][55]. Further protein preparation details can be found in the Supplemental Information [53][54][55][56]. Scheme 2. Synthesis of the thiazolyl-indazole heterocycles IV and V 20 Å. This grid was centered on the previously identified catalytic site of M Pro and used for all subsequent ligand docking using default parameters [12]. Structures of compound I -V were manually built and optimized according to the GAFF force field using Avogadro 1.2 [57,58]. Geometry optimizations in PCM water at the B3LYP/6-31G* level of theory were also performed using Gaussian 16 [59]. Both diastereomers (3aR,6R) and (3aS,6R) of IV were constructed (IV-E and IV-Z respectively). Compound V, with the exocyclic double bond removed, resulted in an altered cyclohexane structure that placed the methyl and isopropyl groups cis (V-E) or trans (V-Z) to each other ( Table 1). As such, V-E and V-Z were also built. Schr€ odinger's Glide Docking program was used to dock these molecules into the catalytic site of M Pro [60,61]. Two poses, the primary pose (1 ) that reports the most favorable GScore and an additional top scoring secondary pose (2 ) with different binding interactions, were selected for further MD evaluation. Docking tools such as Glide are best used as idea generators and selecting poses with different binding interactions helps speed surface coverage in subsequent MD simulation. These poses and their GScores are viewable in Table S1. Additional detail on Glide docking is provided in the Supplemental Information.
As with our treatment of I -V, binding site docking using the Schr€ odinger Glide Docking program was conducted. However, as this program is not intended for peptide docking, this expectedly failed to generate poses aligning with previously described substrate recognition sequencecatalytic site interactions [15,62,64]. Therefore, using UCSF Chimera, each substrate was manually positioned within the M Pro catalytic site to maximize the interactions between the substrate recognition sequence residues P1-Gln and P2-Leu and the M Pro S1 and S2 pockets [65]. The S1 and S2 pockets have been previously noted to be invariably occupied by the P1-Gln and P2-Leu residues, respectively, while low specificity is noted for the P3 and P4 residues [62,64].
Our docking results were used merely as a starting structure for subsequent molecular dynamics simulations. To speed surface coverage for each system studied, we utilized two distinct conformations as simulation starting structures. Ideally, experimental structures (or modifications thereof) are best for initiating MD simulation; however, in cases such as the M Pro :substrate complex, where there is no experimental structure available, subsequent surface sampling imparted by the use of molecular dynamics mitigates concerns associated with our use of manually docked initial structures.
After manual positioning, a minimization calculation was then performed under default settings using UCSF Chimera's Minimize Structure tool. The positions of the P1-Gln and P2-Leu sidechains were fine-tuned using Schr€ odinger's Maestro program [63]. Atoms/sidechains were manually moved, with localized minimizations performed with each movement. The P1-Gln sidechain was positioned to allow for contacts with Phe140, His163, and Gln166 in the S1 pocket, while the P2-Leu sidechain was moved into the S2 pocket. Additionally, the P1-Gln -P1 0 amide bond was positioned between catalytic residues His41 and Cys145. Following these movements, optimization and minimization was performed using Schr€ odinger's Protein Preparation Wizard according to the same steps previously detailed. This resulted in a common binding pose for all three substrates. This pose is detailed in Figure S15B. Additionally, ligand interaction diagrams for each substrate are viewable in Table S2.

Benchmarking: X77 inhibitor preparation and docking
The structure of X77 was extracted from PDB code 6w63. Schr€ odinger's Glide Docking program was used to dock X77 into the catalytic site of M Pro . This was conducted using the same protocol as previously detailed for I -V. Two poses were selected for further MD evaluation. The first pose reported the most favorable GScore and the second pose was obtained from 6w63 by aligning the structure onto chain A of our dimerized 6lu7 model and subsequently removing all 6w63 waters and peptide. These poses are viewable in Table S1.

Molecular dynamics simulations
Unrestrained molecular dynamics (MD) simulations were conducted on M Pro bound to compounds I -V, X77, and three 5-mer substrates, using the AMBER18 suite [66]. The 1 and 2 poses for I -V, and X77, were selected for MD analysis. The manually oriented binding poses detailed above were used to initiate MD for the 5-mer substrate. The ff14SB force field was applied to M Pro and the 5-mer substrates [67]. The program antechamber was used to apply the GAFF force field and AM1-BCC charges to I -V and X77 [57,[68][69][70]. All models were neutralized with Na þ ions and explicitly solvated in a TIP3P unit cell using the program tleap [71]. All simulations were performed using the GPU-accelerated pmemd code of AMBER18 [72,73]. Further details describing the MD protocol can be found in the Supplemental Information. In total, nineteen 1000 ns ensembles were generated (two for each compound I -III, IV-(E,Z), and V-(E,Z), one for each of the three possible 5-mer substrates, and two for X77). The total number of simulations performed is detailed in Table S3.

Molecular dynamics analyses
From each 1000 ns ensemble, all frames in which the ligand either sampled binding positions outside of the catalytic site or dissociated from the protein entirely where removed. All data analysis was conducted on the resultant truncated ensembles including only frames in which the ligand interacts with the catalytic site. As such, percent occurrence and thermodynamic data reported below are relative to each truncated ensemble. This approach allows us to evaluate the specific interaction energies between the active site of the M Pro and the various ligands without the variability associated with normal thermal motion and the random components of the MD algorithm. To avoid any binding overestimation associated with this approach, we also track percent dissociation times for each compound. Trajectory visualization was conducted using UCSF Chimera and UCSF ChimeraX [65,74]. MM-GBSA binding free energy, per-residue decomposition, and normal mode entropy analyses were conducted using the AmberTools MMPBSA.py package [75]. Entropic analysis has been shown to scale binding free energies closer to experimental values while also providing improved comparison of binding affinities across models [76]. This is especially important for comparison between diverse structures like peptide substrates and small molecule inhibitors. Entropy calculations were performed on each truncated ensemble of I and III (with II, IV and V excluded due to experimental and computational shortcomings reported in the Results section) with a 12.5 ns interval. Entropy calculations were also performed on the substrates and the X77 inhibitor. Hydrogen bonding, center-of-mass distance (COM), and root-mean-squared deviation (RMSD) analyses were conducted on I -III using the AmberTools cpptraj module. Potential aromatic -π interactions were screened using literature-based COM distance cutoffs as detailed in Table S4 [77][78][79][80][81][82].
Ensemble averaged binding structures were created for the I -III truncated ensembles using cpptraj RMSD-based clustering. Compounds IV and V were excluded from this analysis due to shortcomings relative to I -III. The RMSD of each frame within each ensemble relative to its respective average binding structure was then calculated. Frames reporting an RMSD value below 1.75 Å were considered to sample the average structure. This cutoff value, intended to account for the dynamic behavior of each structure, was determined from trajectory visualization. Note that the percent occurrence of a structure that undergoes more fluctuation will be underreported. When RMSD analysis alone is deemed insufficient, center-of-mass (COM) analysis is additionally used with cutoff distances selected on a case-by-case basis.
Our simulation protocol is optimized for non-covalent inhibitors and cannot be applied to the covalent inhibitors that have been shown to bind to the SARS-CoV-2 M Pro . As such, a computational comparison to the biochemical assay control compound GC376 could not be conducted. In lieu of this, comparison to X77 and the substrates is provided.

PAINS and QikProp analysis
All compounds passed both PAINS screens with no points of concern. There were no 95% range violations or reactive functional groups reported for I -III, however IV and V generated range violations, and ADME screening indicated a reactive functional group in V (Table 1). Aqueous solubilities for II, III and V-E lie notably close to the logS À6.5 lower limit, while both isomers of IV and V-Z exceed it. Predicted brain/ blood partition coefficients show that I displays the lowest potential for CNS activity while III and IV display the most. All compounds have predicted HERG K þ channel logIC 50 values below À5 and as such are flagged as potential concerns. Compound III and V-Z report the largest potential for HERG K þ channel blockage. It should be noted that this does not necessarily suggest that these or other thiazolyl-indazole compounds are not viable drug candidates. As noted in the QikProp version 6.1 User Manual, failed drug candidates displaying QT-prolongation issues exhibit a large range of IC 50 values, with cisapride (IC 50 6.5 nM, logIC 50 -8.187) and grepafloxacin (IC 50 50000 nM, logIC 50 -4.301) noted as examples. Ultimately, drug compound class, required dosage, and bioavailability determine the allowable IC 50 limit [50]. Taken together, these results suggest that I displays the most favorable drug-like properties.
Property analysis largely eliminates IV-(E,Z) and V-(E,Z) as potential drug candidates, as both display shortcomings relative to I -III. ADME screening indicates that the removal of the hydroxy group from II, to form IV, results in IV having an aqueous solubility outside the 95% range of known drug-like molecules, along with an increased potential for CNS activity and HERG K þ channel blockage. This suggests that the hydroxy group of II may be necessary for maintaining drug-like properties. Likewise, ADME screening indicates that removal of the isopropylcyclohexane double bond in III to produce V results in the cyclohexanelinking imine becoming a potentially reactive carbonyl center (Table 1).

Biochemical evaluation
In order to evaluate the ability of the compounds to disrupt the activity of M Pro , we employed a commercially available kit (see methods). We confirmed the validity of the assay using the known covalent inhibitor GC376 obtaining an IC 50 of 32.5 nM, similar to literature values ( Figure S13) [37]. Since the assay relies on the cleavage of an internally quenched EDANS fluorophore (λ ex 360 nm, λ em 460 nm), we first examined the background fluorescence of the synthesized compounds I, II, IV, and V. As seen in Figure S16, compound IV exhibited high fluorescence when excited at 360 nm so it was excluded from biological testing. Due to solubility issues, we were only able to test compounds II and V up to 100 μM and did not observe inhibition of M Pro activity at these concentrations ( Figure S14). There were no solubility issues with compound I, in agreement with our ADME prediction that I is the most water soluble, while II -V lie close to or even surpass the 95% solubility range of known drug-like molecules, and therefore we were able to characterize the inhibitory activity of I across the full range of solubilities ( Fig. 2, IC 50 92.9 μM).

Molecular dynamics analysis
To better understand the atomistic nature of the experimental results, MD simulations initiated from the structurally distinct 1 and 2 Glide poses (Table S1) of each compound bound to M Pro were conducted, resulting in the production of a 1000 ns ensemble for all of the 1 and 2 Glide poses. Initial data analysis began with the removal of frames where the compound either dissociated from M Pro or sampled a binding position outside of the S1, S1', S2, and/or S4 subsites. A simple analysis based solely on the percentage of frames of I -V that dissociated from M Pro suggest that I, III and V-Z show promise as good binders (Table 2). Increased dissociation of II and IV is likely the result of increased structural rigidity caused by the presence of the indazole group. Trajectory visualization suggests that this increased rigidity results in a diminished ability of II and IV to conform to catalytic site dynamics. To further evaluate binding, we calculated MM-GBSA binding free energies and average structures for I -V using truncated ensembles that included only frames where ligand dissociation did not occur ( Table 2). Taken together, the estimated binding free energies, percent dissociation values, and properties prediction suggest that I, III and V-Z show the most promise as potential inhibitors. This agrees with our experimental results for I; unfortunately, we were unable to synthesize III for experimental testing, and V did not show inhibition at concentrations lower than 100 μM. Limited solubility of V prevented us from testing at higher concentrations. It is noteworthy that our synthesis of V produced both the E and Z isomers and so it is possible that our experimental binding analysis at these concentrations was diminished by the presence of the E isomer. Data for each individual seed and the full non-truncated 1000 ns ensembles is detailed in Tables S5 -S11.
Compound V is not likely to be a viable candidate due to its physical properties; however, a comparison to the results obtained with this molecule provide structural insight into the role of the exocyclic isopropyl moiety on the binding of I and III. While our trajectory analysis of I and III does not provide evidence of any specific interactions between the isopropylcyclohexane double bond and the M Pro catalytic site (see below), the slightly reduced MM-GBSA binding estimations for V suggest that the positioning of the isopropyl group and overall conformation of the cyclohexane ring in I and III affect binding affinity. It is notable that I, with its non-planar isopropyl group, displays a somewhat less favorable binding affinity and higher percent dissociation than III which has a planar isopropyl. Likewise, V-E with a non-planar isopropyl oriented similarly to I, displays a less favorable binding affinity and higher percent dissociation than the V-Z diastereomer with the isopropyl group oriented to the opposite side of the molecule. These results suggest that the orientation of the isopropyl is important and that a cis-oriented isopropyl may be less suitable than a planar or trans-oriented isopropyl.
It is similarly possible that the presence of the double bond in or adjacent to the cyclohexane group affects binding affinity. Geometry optimizations performed on I at the B3LYP/6-31G* level of theory in solvent (PCM water) indicate that this ring is largely planar with a pucker at the isopropyl group. The cyclohexane of III is also largely planar, however the pucker is more pronounced and present at the two cyclohexane carbons between the methyl and isopropyl groups. Keeping in mind that II and V did not produce inhibition in our experimental studies, it is notable that the cyclohexane ring in II is a chair-like whereas in V-E and V-Z the cyclohexane adopts twist boat-and chair-like conformations, respectively (Table 2). Taken together, our experimental binding results and molecular modeling suggests that some planarity in the cyclohexane ring may be important to binding.
Per-residue energy decomposition analysis was conducted on each truncated ensemble of I -III (Table 3). Compound II was included for comparison to I and III, even though II did not show good binding at lower concentrations. The catalytic His41, and Met49, Met165, and Gln189 report highly favorable decomposition energies across all compounds, suggesting their importance for binding interactions. Apart from Gln189, these interactions are largely driven by van der Waals effects.
Visualization of the truncated ensembles suggests close interactions between many of the above residues and the phenyl and thiazole groups of I -III. To further elucidate the importance of this interaction, COM distance analysis was conducted between these aromatic rings and the functional groups of the noted sidechains of I -III (Table S12). This analysis suggests the occurrence of various interactions, most notably π-stacking interactions with His41, and S/π interactions with Met49 and Met165. These compounds display similar per-residue decomposition profiles, and thus the occurrence of similar π interactions is not surprising. Additionally, the potential for OH/π interactions with Ser144, SH/π interactions with Cys145, and anion/π interactions with Glu166 are noted with II and III. These Glu166 interactions resemble Asp anion/π interactions reported in the work of Ellenbarger et al. (face-on packing to the aromatic π-cloud) [83]. The Cys SH/π interaction occurs only with II, while potential Ser OH/π interactions occur only with the phenyl group of II and with the thiazole of III. While the interactions in II suggest promise for the indazole scaffold as the basis of a potential M Pro inhibitor,  Table 3 MM-GBSA per-residue total energy decomposition across all compounds/ensembles. All residues with an average total energy contribution of < À0.5 kcal/ mol are shown. Similarly, if a residue contributed an average total energy of > À0.5 kcal/mol for a particular ensemble, that energy is not reported. The four bolded residues, His41, Met49, Met165, and Gln189 are shown to contribute energy decompositions important for binding across all compounds and truncated ensembles.
Residue Average Total Energy Decomposition (kcal/mol) we do not see experimental inhibition at lower concentrations. In addition to the π/π, S/π, OH/π and anion/π interactions, we also looked at hydrogen bonding patterns for the binding of I. Compound I donates a hydrogen bond to Gln189 for 28.56% of the truncated ensembles, in agreement with the decomposition energies reported in Table 3. Solvent-bridged hydrogen bonds are additionally noted with Glu166 (11.17%) and Gln189 (9.33%). Hydrogen bonding patterns for I -III can be found in the Supplemental Information (Tables S13 -S14).

Average binding structures of I and III
Compound I has experimentally verified M Pro binding behavior suggesting that it may prove useful as an inhibitor scaffold for future development. Compound III is experimentally untested, but our computational results and the similarity of the underlying thiazole backbone structure with I warrants its inclusion for further structural analysis. Within each truncated ensemble of I and III, an average binding structure within the catalytic site was identified through cpptraj analysis. These structures, visualized in Fig. 3 and S17, provide deeper insight to the per-residue decomposition energetics and potential π interactions.
Three average binding structures were identified for I, while two were identified for III. These structures are described as dominant or non-dominant, with the dominant structure reporting the higher percent occurrence. Collectively, the dominant and non-dominant average structures shown in Fig. 3 and S17 represent 93.16 and 60.99% of frames within the I and III ensembles, with the remaining frames sampling disordered binding interactions. RMSD and COM graphs visualizing the percent occurrence of these structures over each ensemble are viewable in Figures S18 -S22. MM-GBSA binding free energy analysis of only the frames sampling these average structures was conducted ( Fig. 3B and Figure S17B). The binding free energy of the I dominant structure closely matches the I average ΔΔG (À26.46 vs. À26.24 kcal/mol). Most notably, both average structures of III report ΔΔG's more favorable than the III average ΔΔG (À30.93 and À29.27, vs. À27.79 kcal/mol). This suggests increased binding affinity of the III average structures over disordered binding interactions.
Differences between the three average structures of I are noted in Fig. 4 and S17. RMSD visualization in Figure S18 indicates that sampling of the I non-dominant structures (22.20 and 14.62%) occurs near the beginning of each simulation, while sampling of the dominant structure occurs later in the trajectory. This transition to the dominant structure occurs even though the non-dominant structure (14.62%) displays an increased binding affinity (À28.46 for non-dominant I vs. À26.46 kcal/ mol for dominant I). The sampling of these structures suggests the possibility of an energetic barrier separating these important conformations. The repeated occurrence of the I dominant structure across individual trajectories suggests that it is entropically favored and explains similarities in per-residue decomposition interaction energies between the two ensembles of I despite structural differences in the poses used to initiate simulation (Table 3). Trajectory visualization indicates an RMSDand COM-indistinguishable variation of the I dominant structure where the (R)-carvone group is oriented inwards towards His41, Met49, and Met165 (Fig. 4B). As this variation orients the phenylthiazole groups away from the three noted residues, His π-stacking and S/π interactions are not possible in this conformation.
It should be noted that in order to generate a classic-textbook, nondynamic pharmacophore, Fig. 3A needs to be combined with 3B, i.e. atomistic physics-based simulations capture the thermal motion and energetic variability associated with real molecules, and as such almost never produce rigid pharmacophores. To obtain such a phamacophore, the conformations in Fig. 3A are combined with the energetics and percent occurrences (or importance weightings) shown in 3B. Such a composite analysis would result in an image for compound I with only one orientation (due to the high percent occurrence (56%) and the overall conformer stability (À26 kcal/mol; second most stable, higher than the lowest energy conformer by only 2 kcal/mol and close enough to be considered approximately isoenergetic). Compound III does in fact sample two distinct orientations as noted by the two conformations with similar weightings and isoenergetics. Pharmacophores with multiple orientations are not unusual [84,85].
Mapping regions of I and III to the M Pro subsites reveals commonalities amongst the average structures of I (Fig. 5A -D). All average structures of I occupy the S2 and S4 subsites. This agrees with our perresidue decomposition (Table 3) and COM distance data (Table S12) suggesting interactions with His41, Met49, Met165, and Gln189. Additionally, the reversed I dominant structure (Fig. 5B), the I non-dominant structure (14.62%) (Fig. 5D), and the III dominant structure (Fig. 5E) project their isopropyl and cyclohexane groups into the S2 pocket. Positioning of the isopropyl group into this pocket is expected due to the S2 subsite favoring occupation by a Leu residue. As discussed above, simulations on V suggest that the positioning of the isopropyl group affects affinity and so we were not surprised to see it appear as a common binding motif in I and III; the ability to occupy the S2 subsite is likely related to the conformational preferences of the various isopropyl moieties. Notably, I is the only compound to significantly interact with the S4 subsite, while the dominant structure of III (Fig. 5E) is the only average structure to bridge both the S1 and S2 subsites. Overall, all average structures demonstrate occupation of the previously noted important subsites for inhibitor binding. Further structural modification and optimization of these scaffolds may allow for additional interactions with various subsites and overall improved binding and inhibition.

Comparison to 5-mer substrates
We performed MD simulations on three 5-mer M Pro substrates to provide a detailed binding comparison with the thiazole compounds. As above, truncated ensembles were produced and MM-GBSA binding free energy and per-residue decomposition analyses were conducted. Binding affinities between the VKLQA and VKLQS substrates appear highly similar, with both reporting average ΔΔG's of À34.88 kcal/mol. The VKLQG substrate reports a notably less favorable binding affinity and a higher percent dissociation from the catalytic site. This is expected, as coronavirus M Pro appears to more readily select for Ser and Ala in the P1 0 position over Gly [86,87]. These energies, reported in Table 4, suggest favorable binding for all three substrates, and on a scale that is relatively comparable to I. (Data for each individual seed is viewable in Table S15.) Normal mode calculations were performed in order to include entropic effects in our binding free energy estimations, and to allow for comparison of I and III to the 5-mer substrates. These results are displayed in Table 5. Expectedly, the 5-mer substrates report greater average entropic penalties than that of the phenylthiazole compounds. Corrected binding free energies, calculated by assuming the MM-GBSA ΔΔG to be roughly equivalent to ΔH and subtracting the entropic penalty, suggests that I and III show greater binding affinities than the 5-mer substrates. Compound III stands out amongst the phenylthiazole compounds, with the primary and secondary truncated ensembles displaying corrected binding free energies of À9.07 and À6.52 kcal/mol, respectively, while compound I displays corrected binding free energies in the range of À5.4 to À5.9 kcal/mol. This further supports the promise of these compounds to serve as scaffolds for further inhibitor development. It should be noted that in these benchmarking tests we are using truncated peptides and as such likely underestimate binding. Also, our simulations do not include any chaperoned folding or other in vivo effects.
Per-residue decomposition and hydrogen bonding analyses reveal similar interactions to those described by Chan et al. [27] Many of the same residues noted to interact with their 11-mer substrates also interact favorably with our 5-mer substrates (Table S16). Interactions with these residues are also noted with I (Table 3). While the VKLQG substrate does not demonstrate any consistent hydrogen bonds, backbonebackbone interactions between Glu166 and P3 Lys are reported with both VKLQA Fig. 4. Compound I dominant structure variations. A. I dominant structure with phenylthiazole oriented inwards towards His41, Met49, and Met165. Potential S/π interactions are displayed, and while a π-stacking interaction with His41 is not shown in this frame, the phenyl group location allows for these interactions across frames sampling this structure. A different orientation of this binding structure is also shown in Fig. 3A. B. I dominant structure with (R)-carvone oriented inwards towards His41, Met49, and Met165. No potential π interactions are noted. C. Overlay of both I average structure variations, with A shown in light green and B shown in dark green. While both of the displayed variations occupy the same space and report low RMSD and COM distance values of 0.86 and 1.92 Å respectively, COM analysis reveals that the phenylthiazole groups are located 6.95 Å apart. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) and VKLQS. This has previously been noted as a primary interaction responsible for holding the substrate in place [27]. Notable interactions are also observed between P2 Leu and His163 as well as P1 Gln and Gly143, Cys145, and His163. Some previously unreported, but less consistent hydrogen bonds are also observed. Most notable of these are interactions between the sidechains of P1 Gln and Ser144 and His41, as well as the P3 Lys backbone and the Gln189 sidechain (Table S17). Peptide:M Pro models are available at https://github.com/Parish-Lab/ sarscov2.

Comparison to X77
Our benchmarking results on the tight-binder X77 suggests that our methodology is appropriate. An average MM-GBSA ΔΔG of À44.16 kcal/ Fig. 5. Compounds I and III average structure interactions with M Pro pocket residues. Yellow indicates S4 pocket residues, cyan indicates S2 pocket residues, red indicates S1 pocket residues, and green indicates S1 0 pocket residues. Each average binding structure is displayed, with the I dominant structure and its reversal shown in A and B, the I non-dominant structures shown in C and D, and the III dominant and non-dominant structures shown in E and F, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)  mol is reported for X77 and the simulation did not produce any dissociated frames indicative of a tight binding ligand. The ensemble initiated from the 6w63 binding pose reports a significantly more favorable MM-GBSA binding energy (-49.53 kcal/mol) than that initiated from the top Glide pose (À38.74 kcal/mol) confirming that experimental structures make the best initial poses for subsequent MD, and that there may be a variety of possible active site orientations with a range of favorability. Binding free energy for each individual seed and per-residue decomposition data is reported in Tables S18 -S19. The entropic contribution and adjusted ΔΔG are reported in Table 5. Our simulations confirm that the enthalpy of X77 binding is more favorable than any of the phenylthiazole scaffolds but less favorable than the various substrates. When entropy is considered, X77 is a markedly more favorable binder than either the substrate or the phenylthiazole compounds, underscoring the importance of rigidity imparted by small-molecule drugs, and the need for further development of the thiazole scaffolds.

Conclusion
Our experimental binding analysis suggests that previously synthesized phenylthiazole compound I is a promising scaffold for future inhibitor development. Computational analysis suggests that I, along with the unsynthesized III, pass PAINS screenings and display ADME properties in-line with known drug-like molecules. Docking analysis predicts these molecules to bind favorably to the SARS-CoV-2 M Pro catalytic site and molecular dynamics simulations suggest prolonged interaction with the enzyme. Subsequent MM-GBSA binding free energy and per-residue decomposition calculations suggest that I and III experience energetically favorable binding driven predominantly by interactions with residues His41, Met49, Met165, and Gln189. Importance is also demonstrated for residues Ser144, Cys145, Glu166, and Asp187 across some of the ensemble-averaged structures. A detailed structural comparison suggests that the orientation of the isopropyl moiety and conformation of the cyclohexane ring may be important for M Pro binding. Trajectory visualization, per-residue decomposition and COM analysis supports the importance of specific π interactions between the phenylthiazole moiety and M Pro . Binding is driven by π-stacking interactions with His41 and possibly S/π interactions with Met49 and Met165. Other notable interactions are supported, including OH/π and SH/π interactions with Ser144 and Cys145 respectively and anion/π interactions with Glu166. An entropically corrected binding comparison to 5-mer substrates with demonstrated experimental affinity further suggests that these compounds show promise as scaffolds for future development as SARS-CoV-2 M Pro inhibitor drugs.