Inhibition of SARS-CoV-2 main protease: A repurposing study that targets the dimer interface of the protein
Coronavirus disease-2019 (COVID-19) was firstly reported in Wuhan, China, towards the end of 2019, and, unfortunately, within a short period of time, emerged as a pandemic. The spread and lethality rates of the COVID-19 have ignited studies that focus on the development of therapeutics for either treatment or prophylaxis purposes.In parallel, drug repurposing studies have also come into prominence. In this study, we aimed at having a holistic understanding of conformational and dynamical changes induced by an experimentally characterized inhibitor on main protease (M-pro) which would enable the discovery of novel inhibitors. To this end, we performed molecular dynamics simulations using crystal structures of apo and α-ketoamide-13b-bound M-pro homodimer. Analysis of trajectories pertaining to apo M-pro revealed a new target site, which is located at the homodimer interface, next to the catalytic dyad. Thereafter, we performed ensemble-based virtual screening by exploiting the ZINC and DrugBank databases and identified three candidate molecules, namely eluxadoline, diosmin, and ZINC02948810 that could invoke local and global conformational rearrangements which were also elicited by α-ketoamide-13b on the catalytic dyad of M-pro. Furthermore, ZINC23881687 was also discerned as a promising candidate due to its interaction with catalytically important residues Glu166 and Ser1. Last but not least, we could find another candidate, namely ZINC20425029, whose mode of action was different. It modulated the dynamical properties of catalytically important residue, Ala285 rather than the catalytic dyad. As such, this study presents valuable findings that might be used in the development of novel therapeutics against SARS-CoV-2 M-pro.
Introduction<!>The simulation protocol<!>Clustering trajectories pertaining to apo M-pro homodimer and binding pocket identification<!>Determination of pharmacophore features and candidate molecule search in the ZINC and DrugBank databases<!>Molecular docking<!>Preparation of ligand-protein complexes and corresponding MD simulation protocols<!>The root-mean-square fluctuation (RMSF) calculations<!>The probability distributions of residue-pair distances<!>Calculation of percent of hydrogen bond occupation<!>Principal component analysis (PCA)<!>Results and discussion<!>The inhibitor induced remarkable conformational rearrangements in the catalytic dyad and substrate-binding pocket of M-pro enzyme<!>Clustering of trajectories pertaining to apo M-pro revealed a possible binding site which is located at the dimer interface
<p>People have been suffering from contagious coronavirus disease-2019 (COVID-19) since December 2019. It is caused by the acute respiratory syndrome coronavirus-2 (SARS-CoV-2) and spreads among people when an infected person makes close contact with uninfected people in the population. [1][2][3][4] Consequently, it has led to a severe pandemic since March 2020.</p><p>Among the most frequently stated symptoms are dry cough, shortness of breath, nasal congestion, abnormal sputum production, and fever as well as debilitation of organs/systems, particularly occurring in the late-stage. 5 Most of the deaths were mainly caused by pneumonia, acute respiratory distress syndrome, and multi-organ failure caused by an exorbitant increase in inflammatory cytokines. 6,7 The extensive genomic studies have shown that SARS-CoV-2 encompasses ca. 30,000 nucleotides and shares 82% sequence identity with the SARS-CoV, hence categorized as a member of Betacoronaviruses [1][2][3]8,9 . Besides high sequence identity, the mechanism of pathogenesis is also similar, where spike protein binds to angiotensin-converting enzyme-2 of the host, which is expressed in the lung, nasopharynx, and intestinal epithelial tissues [10][11][12][13][14] , to enter the cell. After the cellular entry, viral RNA launches translation of large polyproteins, namely polyprotein-1a and -1ab in the ribosome 2,3,[15][16][17] , which are further processed by the M-pro and papain-like protease 18,19 . Moreover, the M-pro of the SARS-CoV-2 shares 96.1% sequence identity with that of the SARS-CoV. It is also important to emphasize that known human proteases possess different cleavage specificity. Therefore, targeting the M-pro enzyme by therapeutics is less likely to pose toxic effects to humans. Consequently, the M-pro has come into prominence as a drug target for preventing viral replication. 1,[20][21][22] In parallel, substantial efforts have been made to unveil 3D structures of the M-pro enzyme 23,24 . It has been revealed that M-pro is a homodimer and each subunit (A and B) is composed of three domains (Domain I, II, and III) along with an N-finger domain that encompasses the first nine residues of the protein. Domain I and II (residues 10 to 99 and 100-182, respectively) adopt a chymotrypsin fold and collectively make up the substrate-binding pocket (SBP). It is located at the interface of the two domains and its conformation is modulated by interactions formed between the N-finger and Glu166 1,20 . The SBP is responsible for the recognition of peptide substrates and their cleavages as well and composed of four subsites, namely S1, S2, S3, and S4 1 . Structural studies have shown that S1 and S2 subsites are well-conserved in Betacoronavirus family [25][26][27][28] . S1 subsite is built up of side-chain atoms of His163 and Phe140 and backbone atoms of Met165, Glu166, and His172 21,29 . It is stabilized by hydrogen bonds formed between residues that are located on subunits A and B of M-pro such as Ser1A-Phe140B, Ser1A-Glu166B, Ser1B-Phe140A, and Ser1B-Glu166A, where A and B denote the name of the subunits 17,30 . Moreover, SBPs of both subunits participate in the formation of the catalytic dyad, which is driven by interactions made between Cys145 and His41  Besides S1, S2 subsite is also crucial as it is responsible for the plasticity of the SBP, which is required for recognition of the protein by its interacting partners 1,32 . Domain III (residues 198 to 303), which is made up of five α-helices, on the other hand, is mainly involved in dimer formation as similar to N-finger domain. 1,[33][34][35] In addition, it also plays a key role in the enzymatic activity of the protein as deletion of the domain caused loss of the enzymatic activity. 1,36 Elucidation of a cysteine residue (Cys145), whose side chain can be readily oxidized, in the catalytic dyad of the M-pro enzyme has prompted studies towards the development of therapeutics that can covalently bind to the residue to inhibit the activity of the enzyme 37,38 .</p><p>To investigate possible binding modes of inhibitors, SARS-CoV-2 M-pro was crystallized with α-ketoamide-11a, 11b, and 13b, as well as antineoplastic carmofur. 1,17,30,39,40 The covalent Michael inhibitors bind to the M-pro by means of nucleophilic attacks made at Cys145. However, this might induce adverse effects, such as allergies, tissue damage, or carcinogenesis, thereby violating their drug safety profiles 41,42 . Consequently, all these attempts have expedited studies that have focused on elaborating the mechanism of Cys145-driven catalytic activity. Kneller et al. has shown that when the distance between Sγ atom of Cys145 and N 2 atom of His41 is 3.8 Å, a hydrogen bond between these atoms would not be formed in the ligand-free SARS-CoV-2 M-pro. In this state of the enzyme, the thiol group of the Cys145 is protonated, whereas the imidazole ring of His41 is neutral. This conceivably indicates that substrate binding or transition state of the enzyme would initiate proton transfer from Cys145 to His41, hence subsequent activation of the catalytic dyad 31 .</p><p>It is also worth mentioning that the M-pro enzyme has been shown to display asymmetry in terms of activation states of the protomers 1 . That is to say, when both protomers of the M-pro are bound with α-ketoamide-13b, only Glu166B-not the one on the subunit A, was shown to be in an inactive state, which was characterized by loss of the H-bond formed between P1 moiety of the α-ketoamide-13b and Glu166B. Moreover, a similar finding was also obtained in a study where SARS-CoV M-pro was shown to display asymmetry upon activation 43,44 .</p><p>The availability of plenty of structural data pertaining to the M-pro enzyme has triggered computer-aided drug design studies for either development of novel therapeutics or repurposing approved molecules. For the latter, virtual screening has been emerged as an effective computational technique for identifying compounds in libraries of small molecules that are most likely to bind to the target protein [45][46][47][48][49][50] .</p><p>In this study, we employed atomistic molecular dynamics (MD) simulations to investigate structural and dynamical properties of apo and α-ketoamide-13b-bound (holo), which is a well-characterized inhibitor, M-pro homodimer. Comparative analyses of the trajectories showed that α-ketoamide-13b increased mobility of the loops found on subunit A and B, on which Glu47 and Arg188 are located, but only disturbing the catalytic dyad of subunit A. On the other hand, it did not induce any change on subunit B. We also clustered the trajectories pertaining to apo M-pro to identify mostly populated conformational states therein. Consequently, we deciphered a possible binding site, which is located at the homodimer interface. The pharmacophore features of this region were modeled and possible candidates, which could bind to the region, were sought in the ZINC 51 and DrugBank [52][53][54][55][56] databases. Thereafter, candidate molecules were docked to the target site on the M-pro and those with favorable binding energies were further tested in molecular dynamics simulations (a total of 16.5 µs) to evaluate their impact on the structure and dynamic of the M-pro.</p><p>Afterward, structural and dynamical properties of ligand/M-pro complexes were compared to those of holo M-pro to identify compounds that could display similar impact as elicited by α-ketoamide-13b. Consequently, diosmin, eluxadoline, and ZINC02948810 were found to be successful candidates. Moreover, ZINC23881687 was considered as a successful candidate even though it did not invoke similar dynamic and structural properties to holo protein, since it formed hydrogen bonds with both Glu166A and Ser1B, which are involved in dimerization. Apart from the aforementioned molecules, we could also identify another candidate, ZINC20425029, which impacted the dynamics of Ala285, which is known to be responsible for the catalytic activity of the enzyme. 1,20 Finally, in addition to local modulation of the Methods Preparation of apo and holo M-pro homodimer systems Crystal structures of apo (PDB ID: 6Y2E) and α-ketoamide-13b-bound M-pro monomer (holo) (PDB ID: 6Y2F) were obtained from the Protein Data Bank (PDB) 1,23,24 . The crystal structure of the apo form of the protein contains all the residues. The homodimer was built up according to the BIOMT information provided in the PDB file using the Visual Molecular Dynamics (VMD) software 57 . Missing residues, Glu47, Asp48, and Gln306, in the crystal structure of holo M-pro, were filled using the GalaxyFill 58 . Afterwards, dimethyl sulfoxide was removed from the crystal, and subsequently, the system was dimerized according to the BIOMT data using the VMD 57 software. The ionization states of residues of apo and holo M-pro dimers were set at pH 7.0 using the ProteinPrep tool of Schrödinger Maestro software 59,60 . Eventually, systems were solvated and neutralized using the Solution Builder of CHARMM-GUI [61][62][63] . The dimensions of the water boxes were set as 104x104x104 Å and 105x105x105 Å for apo and holo systems, respectively. As a side note, the dimensions of the water boxes were determined according to information deposited in the COVID-19 archive of CHARMM-GUI 61 . The systems were neutralized with 0.15 M KCl based on the Monte-Carlo ion displacing method as it is more accurate than distance-based ion displacement 64 . Proteins and water molecules were modeled using the CHARMM36m force-field 65 and TIP3P 66 , respectively.</p><!><p>MD simulations of the apo and holo systems were carried out using the GROningen MAchine for Chemical Simulations (GROMACS) package 67 . The systems were minimized by means of the steepest descent algorithm. Subsequently, the systems were equilibrated in the NVT ensemble for 125 picoseconds (ps), where the Nose-Hoover thermostat was used to maintain the temperature at 310.15 K 68,69 . Furthermore, production runs were done in the NPT ensemble using the Nose-Hoover thermostat and Parrinello-Rahman barostat 70 for a total of 4 µs. In minimization, equilibration, and production steps, the LINCS algorithm was exploited to constraint H-bonds 71 . Two replicates, each of which started with a different velocity distribution, were used for the systems studied.</p><!><p>The distance between Cα atoms of Glu47A and Arg188A was measured in apo M-pro trajectories to determine conformational states pertaining to the SBP. Thereafter, time-line data was converted into the probability plots by setting the frequency range to 1 Å. Consequently, two frames, which represented the most probable states of the SBP, were selected. Subsequently, the selected structures were prepared using the OPLS3e force-field available in the ProteinPrep tool of Schrödinger software 59,60,72 . The possible binding pockets were sought by means of the SiteMap module of Schrödinger, where hydrophobicity, hydrophilicity, druggability, volume, exposure, and enclosure of each possible pocket were calculated and scored, which are collectively given as SiteScore 60,73,74 . The binding sites having close proximity to the homodimer interface and relatively higher SiteScores were selected to be used in drug repurposing studies.</p><!><p>The residues, which are crucial for the function of the M-pro enzyme, were included in pharmacophore modeling. In addition, we also included hotspot residues that were determined using the Hotpoint server 75 . To do so, we compared hotspot residues in both apo and holo trajectories of M-pro homodimer and identified those that displayed lower frequency in holo than in the apo system. The pharmacophore features of each target structure were determined by considering geometrical and chemical properties of the residues Ala116A, Phe140A, Leu141A, Glu166A, His172A, Ser1B, Arg4B, and Gln299B via the Develop Pharmacophore Hypothesis module of Schrödinger 76,77 . Thereafter, the pharmacophore models were assessed according to the relative orientations of the groups included in the model.</p><p>The models which contained groups that were clustered far from each other were excluded to be searched in the ZINC 51 and DrugBank [52][53][54][55][56] databases. Moreover, molecules having less than 700 kDa molecular weight, and 15 rotatable bonds were retrieved from the ZINC and DrugBank databases. Finally, the candidates having a relatively higher match with the pharmacophore model were determined using the Hypothesis Validation tool of Phase 60,78,79 .</p><!><p>The grid files of each target structure were generated to compute interactions between the binding pocket and ligands. To this end, the dimensions of grid boxes were determined according to the centroid of the residues utilized in the pharmacophore model using the Receptor Grid Generation module of Glide 60,[80][81][82] . Before starting molecular docking, the ionization states of ligands were determined at pH 7.0 and stereochemical information was obtained from 3D geometry via the LigPrep module of Schrödinger 60 . Eventually, the ligands were docked to the receptors by using the Glide standard-precision (SP) docking algorithm 60,[80][81][82] . The obtained ligand-protein complexes were assessed according to GScores and orientation of the ligands within the binding pocket. According to that, ligands, which had relatively lower GScores and interacted with at least one of Glu166A, Ser1B, and Arg4B residues, were considered as promising candidates.</p><!><p>The topology and parameter files of the ligands were generated using the Ligand Reader & Modeller of CHARMM-GUI 83 . Subsequently, the obtained topology and parameter files were exploited to prepare the ligand-protein complexes for MD simulations utilizing the Solution Builder of CHARMM-GUI [61][62][63] , where the systems were neutralized with 0.15 M KCl via the Monte Carlo ion displacement. The protein and water molecules were modeled using the CHARMM36m force-field 65 and TIP3P 66 , respectively. The systems were simulated using three replicates per system, each of which started with a different velocity distribution for a total of 16.5 µs following the same procedures explained above.</p><!><p>The RMSF of apo and holo systems were calculated using the 'gmx rmsf' module of GRO-MACS 67 and with following equation:</p><p>, where simulation time over which one wants to average, and coordinates of backbone atom X i at time n were denoted by N, and X i (n), respectively. In this way, the fluctuation pattern of each amino acid was calculated and compared to point out regions where substantial differences were observed.</p><!><p>The distance between certain residues was measured using the 'gmx distance' module of GROMACS 67 . Then, corresponding time-line data was converted into probability plots. To this aim, the minimum and maximum distances were calculated for distance data, and in seriatim, the data was sampled by 2 Å from the minimum distance value to the maximum distance value. Lastly, the frequencies of each interval were calculated and divided by the total frequency, and consequently, the probability of each sampled distance was computed.</p><p>The range of sampling was set to 1 Å only for the selection of frames, which represent different states of the apo system, with higher sensitivity.</p><!><p>In order to calculate the percentage of hydrogen bond occupancy of the ligands with Glu166A and Ser1B, the Hbond plugin of the VMD was utilized 57 . For this, the cut-off distance and angle between acceptor and donor were set to 3.0 Å and 20 degrees, respectively.</p><!><p>In order to explore impacts of candidate molecules on dynamical and structural properties of M-pro homodimer, principal components of the systems, which reflected the collective dynamics of trajectories, were calculated. To do so, the Cα atom of each residue was picked up and the systems were aligned with respect to the Cα atoms of the crystal structure. Subsequently, the eigenvalues and eigenvectors of the systems were calculated by diagonalizing the co-variance matrix of each system using the following equation.</p><p>, where covariance matrix is denoted by C mn as well as M mn ∆r m ∆r n corresponds to a change in position from time-averaged structure for each coordinate of m and n atoms.</p><p>, where set of eigenvalues and eigenvectors of the diagonalized co-variance matrix are shown by δ 2 , and v, respectively.</p><p>Consequently, for calculating 2D PCAs, trajectories of each system were aligned with respect to the first two eigenvalues of holo in 2D space. To this aim, the 'gmx covar' and 'gmx anaeig' modules of the GROMACS were utilized to generate the diagonalized co-variance matrix pertaining to each system and calculate corresponding eigenvalues and eigenvectors 67 . In addition, each system was projected in 2D space with respect to the first eigenvalue to unveil the collective change in each Cα atom throughout the trajectory. To do so, an open-source Python package, namely, ProDy, was used 84,85 .</p><!><p>In this study, we employed atomistic MD simulations using crystal structures of apo and α-ketoamide 13b-bound (holo) SARS-CoV-2 M-pro, where the inhibitor was bound to both subunits A and B of the protein. Trajectories of the systems were compared and conformational changes induced by α-ketoamide 13b were elaborated. Consequently, the most probable conformational states of apo M-pro enzyme were identified and ensemble-based virtual screening was performed on the binding pocket that was identified within the dimer interface. Finally, the capability of candidate molecules to induce conformational changes on the M-pro, which were similar to those induced by α-ketoamide 13b, was interrogated by means of MD simulations.</p><p>The root-mean-square fluctuation (RMSF) profiles of the systems revealed differences in the regions which are known to be responsible for function of SARS-CoV-2 M-pro</p><p>Comparison of RMSF profiles pertaining to the apo and holo systems unveiled differences in the regions that are responsible for the function of the protein. In the presence of αketoamide 13b, fluctuation of N-finger (residues 1 to 9) of subunit A decreased, whereas fluctuation of the short helix on Domain I (residues 54 to 60) and Domain II (residues 100 to 182) were slightly increased (See Figure 1). Herein, it can be thought that modulation of the fluctuation pattern of short helix might also affect the dynamics of the SBP due to their spatial arrangement on the protein. Consequently, this might lead to substantial distinctions in the dynamics of the catalytic dyad of the subunit. In correspondence with that, a recent study has pointed out that residues found in Domain I and II play a pivotal role in modulating the conformation of the SBP as well as the activity of the enzyme 1 .</p><p>As opposed to the above-mentioned dynamical changes in Domain II on subunit A, no change was observed in fluctuation patterns of the same domain of subunit B. On the other hand, the short helix of subunit B displayed increased fluctuation in the holo system. Moreover, the fluctuation pattern of the N-finger on subunit B was not also impacted by the inhibitor, except Ser1B. Since the interaction of Ser1B with Glu166A regulates the conformational state of the SBP 1,20 , we also investigated fluctuation profiles of Glu166A and B, and showed that the former fluctuated more in the presence of the inhibitor while the difference was not pronounced for the latter.</p><!><p>To investigate conformational rearrangements induced by α-ketoamide 13b in the SBP and catalytic dyad, we compared distances between certain atoms: (i) Cα atoms of Glu47 and Arg188 (See Figure 2.c), and (ii) N 2 atom of His41 and Sγ atom of Cys145. The results showed that the inhibitor caused widening of the SBP on subunit A as evidenced by longer distances measured for Glu47-Arg188 pair (See Figure 2.a). In line with this, probability distribution pertaining to the distance measured between His41 and Cys145 showed that the have stemmed from the opening of the loops that bear Glu47 and Arg188 on opposite sites of the SBP (See Figure 3   In addition to the above-mentioned residue-pairs, we also measured the distance between Glu166A and Ser1B since this interaction was shown to impact conformation of the SBP (See Figure 4) 1 . The results showed that the ligand could slightly weaken the interaction between these two atoms as evidenced by longer distances sampled in the presence of the inhibitor (See Figure 4(Compare blue and red lines)).</p><!><p>The trajectories pertaining to apo M-pro were clustered depending on the conformational state of the SBP on subunit A, which was described by the distance between Glu47-Arg188.</p><p>The results showed that the SBP could be found in three conformational states, namely State 1, 2, and 3 (See Figure 5). State 1 and 2 were used in the subsequent steps of the study.</p><p>Figure 5: The probability distribution of the distance which was measured between Cα atom of Glu47 and Arg188 on subunit A.</p><p>First, conformations pertaining to State 1 and 2 were assessed according to their SiteMap scores, and those with close proximity to the dimer interface and higher score were identified to have a possible binding site on the M-pro (See Figure S1). For pharmacophore modeling, we included residues not only those that were in close proximity to the identified pocket and had a central role in the catalytic activity of the enzyme and dimerization such as Glu166A, and Ser1B 1,15,86,87 , but also the ones that showed lower frequency as a hotspot in the presence of the inhibitor, namely, Ala116A, Phe140A, Leu141A, Arg4B, and Gln299B. Moreover, we also included His172 in the pharmacophore model as it forms a hydrogen bond with Glu166, hence, presumably indirectly modulates the conformation of the SBP 1,43 . Therefore, Ala116A, Phe140A, Leu141A, Glu166A, His172A, Ser1B, Arg4B, and Gln299B were included in the pharmacophore models (See Figure S2). Subsequently, pharmacophore features of the binding pockets were evaluated according to the organization of the groups included in the models. If the groups in the models were oriented far from each other it was not considered in the subsequent database search step. Furthermore, the ZINC 51 and DrugBank [52][53][54][55][56] databases were sought to find molecules that could match the properties of the pharmacophore model created. Consequently, a total of 4,174 molecules was retrieved from the databases and docked to the identified pocket on apo M-pro using the Glide SP Docking scheme (See Figure 6.a and 7.a) 60,[80][81][82] . Candidate molecules that could interact with at least one of Glu166A, Ser1B, and Arg4B residues, were considered as promising ones due to substantial roles of these residues in dimerization and enzymatic activity 1,15,[86][87][88][89] . In this regard, six molecules were selected, whose GScores varied between -4.540 to -7.615 (See Table 1 and Figure 6 and 7).   Comparison of trajectories pertaining to ligand/M-pro and α-ketoamide-13b/M-pro complexes reveals similar conformational rearrangements on the substrate-binding pocket and catalytic dyad</p><p>We sought for molecules that could induce similar dynamical and structural rearrangements which were elicited by α-ketoamide-13b by means of atomistic MD simulations. To this end, we calculated and compared probability distribution plots pertaining to distances, which were measured between (i) Cα atoms of Glu47 and Arg188, and (ii) N 2 atom of His41 and Sγ atom of Cys145, and could identify molecules with similar propensities for mediating conformational changes that were induced by α-ketoamide-13b. Here, it is important to emphasize that the successful molecules bound to the novel target site on M-pro without need of any covalent bond formation and influenced the dynamics of the catalytic dyad as opposed to those which were covalently bound to the catalytic dyad 17, 30,91 .</p><p>The analysis of the trajectories showed that eluxadoline, which has been recently approved to be used in the treatment of Irritable Bowel Syndrome with diarrhea 92 , could induce similar conformational rearrangements induced by α-ketoamide-13b. Importantly, it was also proposed as a promising candidate for inhibition of SARS-CoV-2 M-pro in a recent molecular docking study 93 . Specifically, eluxadoline could enable sampling of longer distances by the Glu47-Arg188 pair located on subunit A (See Figure 8  Another compound identified in this study was diosmin, and like eluxadoline, it has been proposed as a candidate molecule in the treatment of SARS-CoV-2 in a couple of recent in silico studies 8,93,94 . It is also important to emphasize that diosmin has been currently tested to be used in a combinatorial therapy with hesperidin for prophylaxis and treatment of COVID-19 95 . Analysis of the trajectories showed that diosmin caused widening of the loops on which Glu47 and Arg188 are located on subunit A. On the other hand, this rearrangement did not cause a significant change on the structure of the catalytic dyad as seen for α-ketoamide 13b on subunit A (See Figure 9.a and c). The expansion of the loops on subunit B caused a slight distortion in the catalytic dyad as opposed to α-ketoamide 13b (See Figure 9.b and d).  As opposed to the above-mentioned ligands, analysis of the trajectories pertaining to ZINC20425029-bound M-pro showed that the ligand did not have a significant impact on the distances measured between Glu47 and Arg188 as well as His41 and Cys145 (See Figure S3). Herein, we thought to further consider this ligand in the candidate short-list as it increased the distance between Ala285 residues (See Figure 11) located on both subunits as a result of occupation of a hydrophobic cavity, which was formed by Met6A, Ala7A, Val125A, Tyr126A, Met6B, Ala7B, Val125B, and Tyr126B, by para-flourophenylmethyl moiety of the ligand. The importance of the Ala285 pair is due to its role in the regulation of catalytic activity of the dimer. Moreover, Ala285A and Ala285B were oriented closer to each other in SARS-CoV-2 than in SARS-CoV due to Thr to Ala mutation 1,20 . Therefore, modulation of the Ala285 pair might prevent formation of a catalytically competent dimer.</p><p>Recalling the GScore values provided in Table 1, that of ZINC23881687 turned out to be smaller than those of other candidates. On the other hand, the ligand was included in the subsequent molecular dynamics studies as its docking pose was found to make contacts with crucial residues such as Glu166, Ser1 and Arg4. Detailed analysis of the trajectories showed that the ligand did not induce significant changes in the distances measured between Glu47 and Arg188 as well as His41 and Cys145 (See Figure S4) similar to ZINC20425029. On the other hand, it could make H-bonds with main and side chains of the Ser1B (See Figure 12).</p><p>Considering the fact that Ser1B is involved in the dimerization of M-pro and modulation of the SBP the perturbation of the residue might help interfere with the function of the M-pro enzyme.</p><p>Finally, it is also crucial to stress that although ZINC39362669 and ZINC02948810 had close GScore values, the former did not have a significant impact on the structural properties discussed for other candidates (See Figure S5).  The global dynamics of M-pro is differently modulated by the ligands and α-ketoamide 13b</p><p>In addition to investigation of local structural/dynamical features of the systems, we also carried out an essential dynamics analysis to explore dominant collective motions of the systems studied. Firstly, we presented 2D projections of the trajectories with respect to the first two eigenvectors of the holo system, which cumulatively accounted for 36% and 33% of the overall motions of subunit A and B, respectively. The results showed that all the ligands, which were found to be successful in terms of local modulation of the M-pro, caused confinement of the conformational space available to subunit A (See Figure 13 and S7). For subunit B, it was shown that all the candidates could confine the spanned space except diosmin and ZINC02948810 (See Figure 13 and S7). Moreover, we also projected trajectories of ligand/M-pro systems on 3D structure of Mpro along with their first eigenvectors and compared them with that of holo (See Figure 14.a) to shed some light on regions that were influenced by the ligand. We observed that α-ketoamide-13b mainly affected dynamics of the loops on which Glu47 and Arg188 were located on both subunits of the protein (See Figure 14.a). In line with this, eluxadoline impacted the same loops (See Figure 14.b) on both subunits. Similary, ZINC23881687 also affected the same loops on both protomers in addition to Domain III of subunit B, which is extensively involved in dimerization to form a functional M-pro homodimer (See Figure 14.c) 1,35,86 . Diosmin substantially altered the collective dynamics of the abovementioned loops located on subunit B as well as Domain III of subunit A (See Figure S8.a). principal component pertaining to ZINC02948810 showed that the ligand caused changes on the loops, on which Glu47 and Arg188 were located, on subunit B, while it did not impact the corresponding loop on subunit A (See Figure S8.b).</p><p>Lastly, ZINC20425029 displayed its effect on the loops found on both protomers as well as Domain IIIs, thus thought of as a promising candidate like other mentioned molecules (See Figure S8.c). In conclusion, analysis of trajectories pertaining to apo M-pro enzyme revealed a novel binding pocket at the dimer interface, which can be targeted by small therapeutic molecules. Indeed, we showed that occupation of the pocket by eluxadoline, diosmin, and ZINC02948810 could alter conformational properties of both SBP and the catalytic dyad similarly as elicited by α-ketoamide-13b. On the other hand, although ZINC20425029 and ZINC23881687 did not impact the SBP and the catalytic dyad, they induced perturbations on regions that were important for the functional activity and dimerization of the enzyme. Specifically, eluxadoline and ZINC23881687 could make stable H-bonds with Glu166A and Ser1B, respectively two of which are involved in the dimerization. Hence, perturbation of inter-molecular hydrogen bonding network of these two residues might help prevent formation of functional M-pro dimers. In parallel to this, it was revealed in global dynamics of ZINC23881687-bound Mpro as higher distortion observed in Domain III of the protein. Interestingly, ZINC20425029 impacted interaction between Ala285 residues located both subunits. Since the Ala285 pair has been shown to be important for dimer formation, perturbation of the interaction in the Ala285 pair might prevent formation of the catalytically competent dimers. Moreover, despite having similar GScore values with ZINC02948810, ZINC39362669 did not impact properties of either the SBP or the catalytic dyad emphasizing the fact that binding score values, which were obtained from molecular docking studies, might not be enough for precise assessment of successful candidates. As opposed to similar local conformational changes induced by ligands and α-ketoamide-13b on the M-pro enzyme, we observed differences in the global dynamics of the protein. Notably, the ligands confined conformational space available to subunit A and B, while α-ketoamide-13b allowed the protein to sample a wider conformational space compared to that of apo M-pro. This might be due to the fact M-pro might be stuck in one of the minima on the energy surface by means of ligands that help maintain M-pro in the catalytically incompetent state.</p>