key: cord-0690806-lwev56my
authors: Pekel, Hanife; Ilter, Metehan; Sensoy, Ozge
title: Inhibition of SARS-CoV-2 main protease: a repurposing study that targets the dimer interface of the protein
date: 2021-04-13
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2021.1910571
sha: 02d8a8be6d1f58e11b6df31906ea5e6d419b6353
doc_id: 690806
cord_uid: lwev56my

Coronavirus disease-2019 (COVID-19) was firstly reported in Wuhan, China, towards the end of 2019, and 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. Herein, 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 stably interacted with catalytically important residues Glu166 and Ser1 and the target site throughout the simulation. However, it gave positive binding energy, presumably, due to displaying higher flexibility that might dominate the entropic term, which is not included in the MM-PBSA method. Finally, ZINC20425029, whose mode of action was different, modulated dynamical properties of catalytically important residue, Ala285. As such, this study presents valuable findings that might be used in the development of novel therapeutics against M(pro.) Communicated by Ramaswamy H. Sarma

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 (Sohrabi et al., 2020; Zhang, Lin, Sun et al., 2020; Zhou et al., 2020) . Consequently, it has led to a severe pandemic since March 2020. 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 (Lai et al., 2020) . 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 (Gao et al., 2021; .

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 (Ngo et al., 2020; Schoeman & Fielding, 2019; Zhang, Lin, Sun et al., 2020) . Besides high sequence identity, the mechanism of pathogenesis is also similar, where spike protein binds to angiotensinconverting enzyme-2 of the host, which is expressed in the lung, nasopharynx, and intestinal epithelial tissues (Hoffmann et al., 2020; Mossel et al., 2008; Uhal et al., 2011; Walls et al., 2020; Xu et al., 2020) , to enter the cell. After the cellular entry, viral RNA launches translation of large polyproteins, namely polyprotein-1a and -1ab in the ribosome (Hilgenfeld, 2014; Morse et al., 2020; Zhou et al., 2020) , which are further processed by the M pro and papain-like protease (Hegyi & Ziebuhr, 2002; Liang et al., 2020) . 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 (Anand et al., 2002; Pillaiyar et al., 2016; Zhang, Lin, Sun et al., 2020) . In parallel, substantial efforts have been made to unveil 3D structures of the M pro enzyme (Berman et al., 2000; Burley et al., 2019) . 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-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 (Anand et al., 2002; Zhang, Lin, Sun et al., 2020) . 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 (Zhang, Lin, Sun et al., 2020) . Structural studies have shown that S1 and S2 subsites are well-conserved in Betacoronavirus family Muramatsu et al., 2016; Su et al., 2020; Yang et al., 2005) . S1 subsite is built up of side-chain atoms of His163 and Phe140 and backbone atoms of Met165, Glu166 and His172 (Li & Kang, 2020; Yang et al., 2003) . 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 . Moreover, the SBPs of both subunits participate in the formation of the catalytic dyad, which is driven by interactions made between Cys145 and His41 (Kneller et al., 2020; Zhang, Lin, Sun et al., 2020) . 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 (Zhang, Lin, Kusov et al., 2020; Zhang, Lin, Sun et al., 2020) .

Domain III (residues 198 to 303), which is made up of five a-helices, on the other hand, is mainly involved in dimer formation as similar to N-finger domain (Huang et al., 2004; Shi & Song, 2006; Zhang, Lin, Sun et al., 2020) . 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 (Shi et al., 2004; Zhang, Lin, Sun et al., 2020) .

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 (Hosseini-Zare et al., 2020; Paasche et al., 2014) . To investigate possible binding modes of inhibitors, SARS-CoV-2 M pro was crystallized with a-ketoamide 11a, 11b and 13b, as well as antineoplastic carmofur Ullrich & Nitsche, 2020; Zhang, Lin, Sun et al., 2020) . 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 (Bauer, 2015; Frecer & Miertus, 2020) . 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 Sc atom of Cys145 and N2 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 (Kneller et al., 2020) .

It is also worth mentioning that the M pro enzyme has been shown to display asymmetry in terms of activation states of the protomers (Zhang, Lin, Sun et al., 2020) . That is to say, when both protomers of the M pro are bound with a-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 a-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 (Chen et al., 2006; Hill & Levitzki, 1980) .

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 (Kapetanovic, 2008; Leelananda & Lindert, 2016; Macalino et al., 2015; Manas & Green, 2017; Melo-Filho et al., 2019; Onawole et al., 2020) .

In this study, we employed atomistic molecular dynamics (MD) simulations to investigate structural and dynamical properties of apo and a-ketoamide 13b-bound (holo), which is a well-characterized inhibitor, M pro homodimer, both of which were used as reference systems in the study. Moreover, we also searched for candidate molecules in ZINC (Koes & Camacho, 2012) and DrugBank (Knox et al., 2011; Law et al., 2014; Wishart et al., 2006 Wishart et al., , 2008 Wishart et al., , 2018 databases which could invoke dynamical and structural changes in the enzyme that was elicited by a-ketoamide 13b. The detailed analyses of trajectories of ligand/M pro complexes showed that five ligands might be used as potential inhibitors to prevent the enzymatic activity.

Crystal structures of apo (PDB ID: 6Y2E) and a-ketoamide 13b-bound M pro monomer (holo) (PDB ID: 6Y2F) were obtained from the Protein Data Bank (PDB) (https://www. rcsb.org/) (Berman et al., 2000; Burley et al., 2019; Zhang, Lin, Sun et al., 2020) . 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 (Humphrey et al., 1996) . Missing residues, Glu47, Asp48, and Gln306, in the crystal structure of holo M pro, were filled using the GalaxyFill (Coutsias et al., 2004) . Afterwards, dimethyl sulfoxide was removed from the crystal, and subsequently, the system was dimerized according to the BIOMT data using the VMD (Humphrey et al., 1996) 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€ odinger Maestro software (Sastry et al., 2013; Schr€ odinger Release, 2018) . Eventually, the systems were solvated and neutralized using the Solution Builder of CHARMM-GUI (Brooks et al., 2009; Jo et al., 2008; Lee et al., 2016) . The dimensions of the water boxes were set to 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 (Jo et al., 2008) . 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 (Wu et al., 2014) . Proteins and water molecules were modeled using the CHARMM36m force-field (Huang et al., 2017) and TIP3P (Mark & Nilsson, 2001) , respectively.

MD simulations of the apo and holo systems were carried out using the GROningen MAchine for Chemical Simulations (GROMACS) package (Abraham et al., 2015) . 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 (Hoover, 1985; Nos e, 1984) . Furthermore, production runs were done in the NPT ensemble using the Nose-Hoover thermostat and Parrinello-Rahman barostat (Parrinello & Rahman, 1981) for a total of 4 ls. In minimization, equilibration, and production steps, the LINCS algorithm was exploited to constraint H-bonds (Hess et al., 1997) . Two replicates, each of which started with a different velocity distribution, were used for the systems studied.

The distance between Ca 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€ odinger software (Roos et al., 2019; Sastry et al., 2013; Schr€ odinger Release, 2018) . The possible binding pockets were sought by means of the SiteMap module of Schr€ odinger, where hydrophobicity, hydrophilicity, druggability, volume, exposure, and enclosure of each possible pocket were calculated and scored, which are collectively given as SiteScore (Halgren, 2007 (Halgren, , 2009 ; Schr€ odinger Release, 2018). The binding sites having close proximity to the homodimer interface and relatively higher SiteScores were selected to be used in drug repurposing studies.

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 (Tuncbag et al., 2010) . 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€ odinger (Loving et al., 2009; Salam et al., 2009) . Thereafter, the pharmacophore models were assessed according to the relative orientations of the groups included in the model. The models which contained groups that were clustered far from each other were excluded to be searched in the ZINC (Koes & Camacho, 2012) and DrugBank (Knox et al., 2011; Law et al., 2014; Wishart et al., 2006 Wishart et al., , 2008 Wishart et al., , 2018 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 (Dixon, Smondyrev, Knoll et al., 2006; ; Schr€ odinger Release, 2018).

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 Schr€ odinger Release, 2018) . 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€ odinger (Schr€ odinger Release, 2018). Eventually, the ligands were docked to the receptors by using the Glide standard-precision (SP) docking algorithm Halgren et al., 2004; Schr€ odinger Release, 2018) . 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.

The topology and parameter files of the ligands were generated using the Ligand Reader & Modeller of CHARMM-GUI (Kim et al., 2017) . 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 (Brooks et al., 2009; Jo et al., 2008; Lee et al., 2016) , 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 (Huang et al., 2017) and TIP3P (Mark & Nilsson, 2001) , 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 ls following the same procedures explained above.

The polar and non-polar contributions to the free energy of the ligand/M pro systems were calculated via the g_mmpbsa module (Baker et al., 2001; Kumari et al., 2014) . To do so, the radii of the ions and temperature were adjusted according to the utilized ions and set temperature in the classical MD simulations, respectively. As a side note, for the calculation of the electrostatic interactions, the dielectric constant was set to 4.

where G solvation, G polar and G non-polar corresponded to the free energy of solvation, electrostatic and non-electrostatic contributing terms to the free energy.

In addition, molecular mechanics energy of the bonded and non-bonded terms were computed as shown in the below.

where non-bonded terms refer to the van der Walls and electrostatic interactions.

Lastly, the free energy of the complex was computed as a sum of the average molecular mechanics potential energy and free energy of solvation, meaning that the entropic contributors were not included in the calculations shown in below.

The root-mean-square fluctuation (RMSF) calculations

The RMSF of apo and holo systems were calculated using the 'gmx rmsf' module of GROMACS (Abraham et al., 2015) and with the following equation:

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.

The RMSD of a-ketoamide 13b, eluxadoline, and ZINC23881687 were calculated by aligning the ligand-bound protein with respect to their first frames in the VMD (Humphrey et al., 1996) . Thereafter, the RMSD profiles of the above-mentioned systems were computed according to the following formula:

where the number of total atoms and the position of atom n at time t were denoted by N and d n (t), respectively.

The distance between certain residues was measured using the 'gmx distance' module of GROMACS (Abraham et al., 2015) . 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. 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.

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 (Humphrey et al., 1996) . For this, the cut-off distance and angle between acceptor and donor were set to 3.0 Å and 20 degrees, respectively.

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 Ca atom of each residue was picked up and the systems were aligned with respect to the Ca 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.

where covariance matrix is denoted by C mn as well as M mn Dr m Dr n corresponds to a change in position from time-averaged structure for each coordinate of m and n atoms.

where set of eigenvalues and eigenvectors of the diagonalized co-variance matrix are shown by d 2 , and v, respectively. 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 (Abraham et al., 2015) . In addition, each system was projected in 2D space with respect to the first eigenvalue to unveil the collective change in each Ca atom throughout the trajectory. To do so, an open-source Python package, namely ProDy, was used (Bakan et al., 2011 (Bakan et al., , 2014 .

In this study, we employed atomistic MD simulations using crystal structures of apo and a-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 a-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 a-ketoamide 13b, was interrogated by means of MD simulations.

The root-mean-square fluctuation (RMSF) profiles of the systems revealed differences in the regions which are known to be responsible for the function of SARS-CoV-2 M pro 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 a-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 (Zhang, Lin, Sun et al., 2020) .

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 (Anand et al., 2002; Zhang, Lin, Sun et al., 2020) , 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.

To investigate conformational rearrangements induced by a-ketoamide 13b in the SBP and catalytic dyad, we compared distances between certain atoms: (i) Ca atoms of Glu47 and Arg188 (See Figure 2c) , and (ii) N2 atom of His41 and Sc 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 2a ). In line with this, probability distribution pertaining to the distance measured between His41 and Cys145 showed that the ligand caused a loss in the coordination of the catalytic dyad on subunit A, which might have stemmed from the opening of the loops that bear Glu47 and Arg188 on opposite sites of the SBP (See Figure 3a) . On the other hand, the inhibitor led to neither distortion of the SBP nor coordination loss in the catalytic dyad on subunit B (See Figure 3b) . Interestingly, Glu47 and Arg188 could come closer in the presence of the inhibitor on subunit B (See Figure 2b) .

Herein, it is important to emphasize that the same asymmetrical behavior shown for the protomers of the Figure 6 . (a) The identified binding pocket was shown on 3D structure of M pro in surface representation, which adopted State 1. The protein was shown in New Cartoon representation. The binding poses of (b) ZINC02948810, (c) ZINC39362669, (d) diosmin, and (e) eluxadoline within the pocket were also shown (Stone, 1998). a-ketoamide 13b-bound enzyme was also observed in energy values calculated between the inhibitor and the protomers of the enzyme (See Methods for details). While the inhibitor was stably bound to one of the subunits during the course of the simulation, the other one left the active site in both replicates of the system. In the first trajectory, the ligand left the subunit B around 63 ns while it detached from the same subunit around 200 ns in the second trajectory (See Supplementary material Figure S1 ). Overall, these findings are in line with the observation that protomers act as an asymmetrical activated dimer where only one subunit resembles the properties of the active enzyme.

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) (Zhang, Lin, Sun et al., 2020) . 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)).

Clustering of trajectories pertaining to apo M pro revealed a possible binding site which is located at the dimer interface

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. 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.

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 Supplementary material Figure S2 ). 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 (Stone, 1998) . 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 (Chen et al., 2006; Zhang, Lin, Sun et al., 2020) . Therefore, Ala116A, Phe140A, Leu141A, Glu166A, His172A, Ser1B, Arg4B and Gln299B were included in the pharmacophore models (See Supplementary material Figure  S3 ). 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. 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 Figures 6a and 7a ) Schr€ odinger Release, 2018) . 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 (Douangamath et al., 2020; Goyal & Goyal, 2020; Hilgenfeld, 2014; Zhang, Lin, Sun et al., 2020) . In this regard, six molecules were selected, whose GScores varied between À4.540 to À7.615 (See Table 1 and Figures 6 and 7) .

Comparison of trajectories pertaining to ligand/m pro and a-ketoamide 13b/M pro complexes reveals similar conformational rearrangements on the substratebinding pocket and catalytic dyad

We sought for molecules that could induce similar dynamical and structural rearrangements which were elicited by a-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) Ca atoms of Glu47 and Arg188, and (ii) N2 atom of His41 and Sc atom of Cys145, and could identify molecules with similar propensities for mediating conformational changes that were induced by a-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 Vuong et al., 2020) . 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 (Simr en & Figure 10 . The probability distribution pertaining to distance which was measured between Ca atoms of Glu47 and Arg188 on subunit (a) A and (b) B. The probability distribution pertaining to distance which was measured between N2 atom of His41 and Sc atom of Cys145 on subunit (c) A and (d) B. Tack, 2018), could induce similar conformational rearrangements induced by a-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 . Specifically, eluxadoline could enable sampling of longer distances by the Glu47-Arg188 pair located on subunit A (See Figure 8a) , thereby facilitating distortion of the catalytic dyad therein, which was evidenced by relatively longer distances sampled by His41-Cys145 pair (See Figure 8c ). Interestingly; however, although it did not impact dynamics of the Glu47-Arg188 pair on subunit B (See Figure 8b) , it could slightly distort the catalytic dyad which was evidenced by right shift in the distance probability distribution (See Figure 8d ).

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 (Arun et al., 2020; Chen et al., 2020; Ngo et al., 2020) . 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 (Haggag et al., 2020) . 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 a-ketoamide 13b on subunit A (See Figure 9a and c). The expansion of the loops on subunit B caused a slight distortion in the catalytic dyad as opposed to a-ketoamide 13b (See Figure 9b and d).

Moreover, ZINC029498810 also led to an increase in the distances measured between Glu47 and Arg188 on both subunits (See Figure 10a and b) ; however, this conformational rearrangement did not distort the catalytic dyad on subunit A (See Figure 10c) , whereas it remarkably distorted the catalytic dyad on subunit B (See Figure 10d) .

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 Supplementary material Figure S4 ). Herein, we thought to further consider this ligand in the candidate short-list as, unlike other studied ligands, it increased the distance between Ala285 residues (See Figure 11 and Supplementary material Figure S7 ) 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 (Anand et al., 2002; Zhang, Lin, Sun et al., 2020) . Therefore, modulation of the Ala285 pair might prevent formation of a catalytically competent dimer.

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 Supplementary material Figure S5 ) similar to ZINC20425029. On the other hand, it could make H-bonds with main and side chains of the Ser1B (See Figure  12) . 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.

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 Supplementary material Figure S6 ).

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 Supplementary material Figure S8 ). For subunit B, it was shown that all the candidates could confine the spanned space except diosmin and ZINC02948810 (See Figure 13 and Supplementary material Figure S8 ).

Moreover, we also projected trajectories of ligand/M pro systems on 3D structure of M pro along with their first eigenvectors and compared them with that of holo (See Figure  14a ) to shed some light on regions that were influenced by the ligand. We observed that a-ketoamide 13b mainly affected dynamics of the loops on which Glu47 and Arg188 were located on both subunits of the protein (See Figure  14a ). In line with this, eluxadoline impacted the same loops (See Figure 14b ) 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 14c ) (Douangamath et al., 2020; Shi & Song, 2006; Zhang, Lin, Sun et al., 2020) . Diosmin substantially altered the collective dynamics of the above-mentioned loops located on subunit B as well as Domain III of subunit A (See Supplementary material Figure S9a ). As it is consistent with the mentioned local alterations on the catalytic dyads, the first 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 Supplementary material Figure S9b) . 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 Supplementary material Figure S9c ).

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 a-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 M pro 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 a-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 a-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. Last but not least, all the ligands, except ZINC23881687, gave negative binding energies with M pro meaning that the drug candidates could favorably interact with the target binding site during the course of the simulation (See Supplementary material Table S1 ). On the other hand, ZINC23881687 gave positive binding energy presumably having relatively higher conformational flexibility compared to other ligands studied (See Supplementary material Figure S10 ).

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers

Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra alpha-helical domain

Coronavirus main proteinase (3CLpro) structure: Basis for design of anti-SARS drugs

Drug repurposing against SARS-CoV-2 using E-pharmacophore based virtual screening, molecular docking and molecular dynamics with main protease as the target

Evol and ProDy for bridging protein sequence evolution and structural dynamics

ProDy: Protein dynamics inferred from theory and experiments

Electrostatics of nanosystems: Application to microtubules and the ribosome

Covalent inhibitors in drug discovery: From accidental discoveries to avoided liabilities and designed therapies

The Protein Data Bank

CHARMM: The biomolecular simulation program

RCSB Protein Data Bank: Biological macromolecular structures enabling research

Only one protomer is active in the dimer of SARS 3C-like proteinase

Prediction of the SARS-CoV-2 (2019-nCoV) 3C-like protease (3CLpro) structure: Virtual screening reveals velpatasvir, ledipasvir, and other drug repurposing candidates

Mutation of Glu-166 blocks the substrate-induced dimerization of SARS coronavirus main protease

Quaternary structure of the severe acute respiratory syndrome (SARS) coronavirus main protease

A kinematic view of loop closure

Structure-based design of antiviral drug candidates targeting the SARS-CoV-2 main protease

PHASE: A novel approach to pharmacophore modeling and 3D database searching

PHASE: A new engine for pharmacophore perception, 3D QSAR model development, and 3D database screening: 1. Methodology and preliminary results

Crystallographic and electrophilic fragment screening of the SARS-CoV-2 main protease

Antiviral agents against COVID-19: Structure-based design of specific peptidomimetic inhibitors of SARS-CoV-2 main protease

Glide: A new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy

Extra precision glide: Docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes

Cytokine storm syndrome in coronavirus disease 2019: A narrative review

Targeting the dimerization of the main protease of coronaviruses: A potential broad-spectrum therapeutic strategy

Is hesperidin essential for prophylaxis and treatment of COVID-19 Infection? Medical Hypotheses

New method for fast and accurate binding-site identification and analysis

Identifying and characterizing binding sites and assessing druggability

Glide: A new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening

Potential of coronavirus 3C-like protease inhibitors for the development of new anti-SARS-CoV-2 drugs: Insights from structures of protease and inhibitors

Conservation of substrate specificities among coronavirus main proteases

LINCS: A linear constraint solver for molecular simulations

From SARS to MERS: Crystallographic studies on coronaviral proteases enable antiviral drug design

Subunit neighbor interactions in enzyme kinetics: Half-of-the-sites reactivity in a dimer

SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor

Canonical dynamics: Equilibrium phase-space distributions

Targeting severe acute respiratory syndrome-coronavirus (SARS-CoV-1) with structurally diverse inhibitors: A comprehensive review

3C-like proteinase from SARS coronavirus catalyzes substrate hydrolysis by a general base mechanism

CHARMM36m: An improved force field for folded and intrinsically disordered proteins

VMD: Visual molecular dynamics

Structure of M pro from SARS-CoV-2 and discovery of its inhibitors

Structural basis for the inhibition of SARS-CoV-2 main protease by antineoplastic drug carmofur

CHARMM-GUI: A web-based graphical user interface for CHARMM

Computer-aided drug discovery and development (CADDD): In silico-chemico-biological approach

CHARMM-GUI ligand reader and modeler for CHARMM force field generation of small molecules

Structural plasticity of SARS-CoV-2 3CL M pro active site cavity revealed by room temperature X-ray crystallography

DrugBank 3.0: A comprehensive resource for 'omics' research on drugs

ZINCPharmer: Pharmacophore search of the ZINC database

g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

Extrarespiratory manifestations of COVID-19

DrugBank 4.0: Shedding new light on drug metabolism

CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field

Computational methods in drug discovery

Progress in developing inhibitors of SARS-CoV-2 3C-like protease. Microorganisms

Interaction of the prototypical a-ketoamide inhibitor with the SARS-CoV-2 main protease active site in silico: Molecular dynamic simulations highlight the stability of the ligandprotein complex

Energetic analysis of fragment docking and application to structure-based pharmacophore hypothesis generation

Role of computeraided drug design in modern drug discovery

CADD medicine: Design is the potion that can cure my disease

Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K

Discovery of new potent hits against intracellular Trypanosoma cruzi by QSAR-based virtual screening

Learning from the Past: Possible Urgent Prevention and Treatment Options for Severe Acute Respiratory Infections Caused by 2019-nCoV

SARS-CoV replicates in primary human alveolar type II cell cultures but not in type I-like cells

SARS-CoV 3CL protease cleaves its C-terminal autoprocessing site by novel subsite cooperativity

Computational Determination of Potential Inhibitors of SARS-CoV-2 Main Protease

A unified formulation of the constant temperature molecular dynamics methods

COVID-19: CADD to the rescue

Evidence for substrate binding-induced zwitterion formation in the catalytic Cys-His dyad of the SARS-CoV main protease

Polymorphic transitions in single crystals: A new molecular dynamics method

An overview of severe acute respiratory syndrome-coronavirus (SARS-CoV) 3CL protease inhibitors: Peptidomimetics and small molecule chemotherapy

OPLS3e: Extending force field coverage for drug-like small molecules

Novel method for generating structure-based pharmacophores using energetic analysis

Protein and ligand preparation: Parameters, protocols, and influence on virtual screening enrichments

Coronavirus envelope protein: Current knowledge

4: Glide. Schr€ odinger, LLC

The catalysis of the SARS 3C-like protease is under extensive regulation by its extra domain

Dissection study on the severe acute respiratory syndrome 3C-like protease reveals the critical role of the extra domain in dimerization of the enzyme: Defining the extra domain as a new target for design of highly specific protease inhibitors

New treatments and therapeutic targets for IBS and other functional bowel disorders

World Health Organization declares global emergency: A review of the 2019 novel coronavirus (COVID-19)

An efficient library for parallel ray tracing and animation

CoV-2 activities in vitro of Shuanghuanglian preparations and bioactive ingredients

HotPoint: Hot spot prediction server for protein interfaces

Regulation of alveolar epithelial cell survival by the ACE-2/angiotensin 1-7/ Mas axis

The SARS-CoV-2 main protease as drug target

Feline coronavirus drug inhibits the main protease of SARS-CoV-2 and blocks virus replication

Structure, Function, and Antigenicity of the SARS-CoV

DrugBank 5.0: A major update to the DrugBank database

DrugBank: A knowledgebase for drugs, drug actions and drug targets

DrugBank: A comprehensive resource for in silico drug discovery and exploration

Risk factors associated with acute respiratory distress syndrome and death in patients with coronavirus disease

CHARMM-GUI Membrane Builder toward realistic biological membrane simulations

A new coronavirus associated with human respiratory disease in China

Biological efficacy of perpendicular type-I collagen protruded from TiO2-nanotubes

Design of wide-spectrum inhibitors targeting coronavirus main proteases

The crystal structures of severe acute respiratory syndrome virus main protease and its complex with an inhibitor

Ketoamides as Broad-Spectrum Inhibitors of Coronavirus and Enterovirus Replication: Structure-Based Design, Synthesis, and Activity Assessment

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved a-ketoamide inhibitors

A pneumonia outbreak associated with a new coronavirus of probable bat origin

The authors thank Istanbul Medipol University for providing access to the High-Performance Computing System so as to run some of the MD simulations. The rest of the numerical calculations reported in this paper were performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).

The authors declare that there is no conflict of interest.

This study has not been supported by any funding agencies/companies.