Lead Discovery of SARS-CoV-2 Main Protease Inhibitors through Covalent Docking-Based Virtual Screening.

During almost all 2020, coronavirus disease 2019 (COVID-19) pandemic has constituted the major risk for the worldwide health and economy, propelling unprecedented efforts to discover drugs for its prevention and cure. At the end of the year, these efforts have culminated with the approval of vaccines by the American Food and Drug Administration (FDA) and the European Medicines Agency (EMA) giving new hope for the future. On the other hand, clinical data underscore the urgent need for effective drugs to treat COVID-19 patients. In this work, we embarked on a virtual screening campaign against the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) Mpro chymotrypsin-like cysteine protease employing our in-house database of peptide and non-peptide ligands characterized by different types of warheads acting as Michael acceptors. To this end, we employed the AutoDock4 docking software customized to predict the formation of a covalent adduct with the target protein. In vitro verification of the inhibition properties of the most promising candidates allowed us to identify two new lead inhibitors that will deserve further optimization. From the computational point of view, this work demonstrates the predictive power of AutoDock4 and suggests its application for the in silico screening of large chemical libraries of potential covalent binders against the SARS-CoV-2 Mpro enzyme.

Introduction

At the end of 2019, a new virus belonging to the coronaviridae family initiated an epidemic of pulmonary disease in Wuhan, the capital of the Hubei province in China, and has since spread worldwide.[1] The new coronavirus has been called severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), considering its similarity to the first SARS coronavirus (SARS-CoV), while the disease has been named coronavirus disease 2019 (COVID-19). On March 11th, 2020, the World Health Organization (WHO) declared the pandemic outbreak. As of December 28th, 2020, the WHO reports more than 75 million reported cases and more than 1.6 million deaths worldwide.[2] The Americas and Europe are the regions most affected by the pandemic in terms of the number of confirmed cases and deaths, while Africa and Western Pacific areas the least affected. Moreover, all over the Americas, most of Europe, and Africa, the transmission is classified as “community transmissions” which means that these regions are now experiencing large outbreaks of local transmission.[2] At the beginning of the COVID-19 outbreak, most of the cases were connected to infections contracted after eating animals from the Wuhan market. In a short time, the virus’s high efficiency in spreading among people led to an exponential growth rate, culminating in the COVID-19 pandemic. In this scenario, the dire need for vaccines or for specific agents helpful in curing COVID-19 and controlling its symptoms induced the research community worldwide to put considerable efforts into the identification of possible SARS-CoV-2 druggable targets and in discovering the first agents to modulate them.[3] Indeed, at the end of 2020, these efforts have culminated with the approval of vaccines by the American Food and Drug Administration (FDA) and the European Medicines Agency (EMA) giving new hope for the future. However, the inability to vaccinate the world’s population in a short amount of time and, at the same time, the emergence of new variants of the wild-type strain of SARS-CoV-2 that could hamper the efficacy of the developed vaccines put in evidence the urgent need to develop the first effective drugs to treat COVID-19 patients.

The SARS-CoV-2 genome consists of about 30,000 nucleotides encoding for two polyproteins, namely, pp1a and pp1ab, which are proteolytically cleaved by the main protease (Mpro) and papain-like protease (PLpro) in both structural and non-structural viral proteins.[4]

SARS-CoV-2 Mpro is a chymotrypsin-like cysteine protease consisting of two monomers, each of which includes three domains (I, II, and III), with the catalytic dyad C145–H41 located between the domains I and II.[5] The Mpro crystal structure in complex with the peptide inhibitor N3 was first deposited in the Protein Data Bank (PDB code: 6LU7).[5] Since this first deposition, more than 190 Mpro structures were deposited in the PDB in complex with several peptides, low-molecular-weight molecules as well as fragments. Similar to PLpro,[6] the interest in targeting this enzyme is fueled by its critical role in the maturation of key viral enzymes, which explains why Mpro inhibition blocks viral replication. This enzyme targets at least 11 sites on the 1ab polyprotein and recognizes the Leu-Gln-(Ser, Ala, Gly) sequences, cleaving the peptide bond between Gln and the adjacent amino acid. Since no human proteases are known to recognize this sequence, the selective inhibition of Mpro should avoid off-target effects and, indeed, a series of reversible and irreversible inhibitors have already been described.[5−14] Moreover, several theoretical approaches (i.e., structure- and ligand-based virtual screening, free-energy, and molecular dynamics calculations) have been demonstrated to be useful in the identification of Mpro negative modulators.[15−23]

Our group has a consolidated experience in the design, synthesis, and biological evaluation of pseudopeptides or peptidomimetics acting as cysteine protease inhibitors as valid therapeutic agents for the treatment of infectious diseases as well as in the application of receptor-based virtual screening for the discovery of biologically active compounds. Thus, in the present work, we decided to give our contribution to the field by undertaking a receptor-based VS campaign of our in-house focused chemical library of cysteine protease inhibitors to cherry-pick the most promising candidates to submit to the enzyme inhibition assays. Our starting chemical database results from our longstanding efforts to develop small pseudopeptides and conformationally constrained peptidomimetics able to inactivate several proteases and characterized by different types of warheads acting as Michael acceptors, such as vinyl sulfones, vinyl esters, vinyl amides, vinyl ketones, vinyl phosphonates, and vinyl nitriles, all able to form covalent adducts with the active site thiol function.[24−33] Given the nature of this focused library, we decided to use a covalent docking protocol employing the so-called “flexible side-chain method” available within the AutoDock4 (AD4) docking software.[34] In this method, the ligand is attached in a random conformation to the target residue allowing to model the covalently bound ligand as a flexible side chain in the AD4 simulation. The choice of this protocol was dictated by the high accuracy in reproducing the native binding conformations of covalent ligands’ experimental poses in benchmarking docking experiments. Moreover, the same method has already been successfully used by us in previous studies on some of the ligands present in our in-house focused chemical library against the cathepsin L-like enzyme rhodesain of Trypanosoma brucei rhodesiense.[24−26]

Through these VS experiments, we selected 15 compounds for their potential binding against the SARS-CoV-2 Mpro, which were validated in vitro for their inhibitory activity against the enzyme. This strategy allowed us to identify two new lead compounds, with the most active ones being further profiled for their ability to act as an irreversible Michael acceptor.

Computational and Experimental Methods

Covalent Docking VS

AD4[35,36] was employed for the docking calculations. To probe the formation of the covalent adduct, a specific docking protocol devised by Bianco et al., namely, the “flexible side-chain method”, was employed.[34] This protocol needs to adapt the residue taking part in the covalent bond by attaching the ligand to its side chain; this modified residue is then considered flexible during the docking calculation. To this end, using the Maestro suite,[37] we modeled all the ligands present in the in-house focused library with two extra atoms where the nucleophilic attack by the reactive thiol would take place, specifically, a sulfur and a carbon atom, to match the corresponding atoms in C145. The crystal structure of SARS-CoV-2 Mpro from the Protein Data Bank (PDB code 7BQY)[5] was downloaded and prepared for docking using the Protein Preparation Wizard tool within Maestro. Then, with the aid of the scripts provided by the AD4 website,[38] the ligands (106 compounds from our in-house database, described in the Results and Discussion section and provided in SI, as well as boceprevir and telaprevir)[39] and the C145 residue were overlapped. Subsequently, the receptor grid maps were calculated with the AutoGrid4 software, mapping the receptor interaction energies using the ligand atom types as probes. The grid of 60 Å × 60 Å × 60 Å with 0.375 Å spacing was centered on the coordinates of the ligand originally present in the 7BQY crystal.[5] Finally, the actual docking was attained for each ligand separately, keeping the remodeled C145 residue as flexible. This permitted to sample the torsional flexibility of each ligand within the Mpro enzyme. For the docking simulations, the Lamarckian Genetic Algorithm (LGA) was used. Given the high torsional flexibility, 200 independent docking runs of LGA were attained for each compound. Each docking run consisted of 20 million energy evaluations using the Lamarckian genetic algorithm local search (GALS) method. The GALS method assesses a population of possible docking poses and propagates the most successful entities from each generation into the subsequent one. A low-frequency local search according to the method of Solis and Wets was applied to docking runs to ensure that the resulting solution represents a local minimum. All dockings described in this work were performed selecting a population size of 150, and 300 rounds of Solis and Wets local search were applied with a probability of 0.06. A mutation rate of 0.02 and a crossover rate of 0.8 were set to generate new docking trials for subsequent generations, and the best individual from each generation was propagated over the next generation. All the other settings were left at their default value. The docking results from each calculation were clustered based on root-mean-square deviation (rmsd) (solutions differing by less than 2.0 Å) between the cartesian coordinates of the atoms and were ranked based on the predicted binding free energy (ΔGAD4). All the images were rendered using the UCSF ChimeraX Molecular Modeling Software.[40]

Molecular Dynamics Simulations

All-atom molecular dynamics (MD) simulations were performed on the 6/Mpro and 10/Mpro complexes (for the structures of compounds 6 and 10 see Table of the Results and Discussion section) obtained by the AD4 software using the Desmond module[41,42] of the Schrödinger software package. The system builder panel was used to set the initial system for MD calculation. The complexes were embedded in a parallelepiped box and solvated with TIP3P water models,[43] and the −3 negative charge of the system was neutralized using 3 Na+ ions. The systems were equilibrated employing the NPT ensemble using the default Desmond protocol, and it included a total of eight steps, among which the first 7 were short simulations (equilibration phase) steps at increasing temperature and decreasing restraints on the solute. The equilibrated systems were then subjected to a 100 ns MD production run with PBC conditions and NPT ensemble using the OPLSe forcefield.[44] During the simulation, 1 atm pressure and 300 K temperature of the system were set employing a Martyna–Tobias–Klein barostat[45] and Nose–Hoover chain thermostat.[46]

Table 1

Structures of the Selected Dipeptidyl Derivatives along with the Predicted ΔGAD4 of the Complex Formation

Protein Expression and Purification

The pMal plasmid harboring the C-terminal hexahistidine-tagged sequence of SARS-CoV-2 Mpro was kindly provided by Prof. John Ziebuhr (Justus Liebig University Gießen, Germany). The sequence contained the native nsp4/nsp5 Mpro cleavage site between MBP and Mpro as well as the native nsp5/nsp6 cleavage site between Mpro and the hexahistidine tag, thus enabling the purification of native Mpro. Protein expression was carried out in Escherichia coli strain BL21-Gold (DE3) (Agilent Technologies, Santa Clara, CA, USA). Cells were grown in LB medium with the corresponding antibiotic (100 μg/mL ampicillin) at 37 °C to an OD600 of ∼0.5 and induced with 0.3 mM isopropyl-d-thiogalactoside. Proteins were expressed at 18 °C for 16 h and harvested by centrifugation. For purification, cells were resuspended in lysis buffer (20 mM Tris–HCl pH 7.8, 150 mM NaCl, 20 mM imidazole) and lysed by sonication (Sonoplus; Bandelin, Berlin, Germany). The lysate was cleared by centrifugation, and the supernatant was immediately subjected to immobilized metal affinity chromatography (IMAC) on a HisTrap HP 5 ml column. After washing with 20 column volumes of buffer A (20 mM Tris–HCl pH 7.8, 200 mM NaCl, 20 mM imidazole), the protein was eluted in a linear gradient of buffer B (20 mM Tris–HCl pH 7.8, 200 mM NaCl, 500 mM imidazole). To suppress early cleavage of the hexahistidine tag, IMAC was performed with buffers cooled on ice. After IMAC, the eluted protein was subjected to a gel filtration step (HiLoad 16/600 Superdex 75 pg column; GE Healthcare, Chicago, IL, USA) in SEC buffer [20 mM Tris–HCl pH 7.8, 150 mM NaCl, 1 mM ethylenediaminetetraacetic acid (EDTA), 1 mM dithiothreitol (DTT)]. Prior to shock freezing in liquid nitrogen, the eluted protein was diluted to 10 μM and 10% (v/v) glycerol was added. Throughout purification, protein concentrations were measured via absorbance at 280 nm using a NanoDrop 2000 Spectrophotometer (Thermo Scientific Waltham, Massachusetts). Sample purity was assessed via Coomassie brilliant blue-stained sodium dodecyl sulfate–polyacrylamide gel electrophoresis.

Fluorometric Enzyme-Activity Studies

Fluorometric enzyme assays were performed on a TECAN Spark 10M (Agilent Technologies, Santa Clara, USA) using Dabcyl-KTSAVLQSGFRKME-Edans (Genescript, New Jersey, USA) as a FRET-substrate at 5 μM. Released Edans was excited at a wavelength of 335 nm (slit 20 nm), and fluorescence was recorded at 493 nm (slit 20 nm). The assay was carried out with 250 nM of isolated Mpro at 25 °C in a final volume of 200 μL: to 170 μL buffer (20 mM Tris–HCl buffer pH 7.5, 200 mM NaCl 0.1 mM EDTA, 1 mM DTT), 20 μL of inhibitor solutions in dimethyl sulfoxide (DMSO) or pure DMSO as control and 5 μL of a 10 μM Mpro solution in SEC buffer were added. The reaction was initiated by the addition of 5 μL of 200 nM substrate in DMSO under rigorous mixing and monitored for 10 min if not stated otherwise. IC50 values as well as kinetic parameters were calculated using GRAFIT (Version 5.0.13; Erithacus Software Limited, East Grinstead, West Sussex, UK). KM determination was performed with varying substrate concentrations and addition of DMSO instead of inhibitor solutions.

Dialysis Assay

Dialysis experiments with SARS-CoV-2 Mpro were performed using a custom-built dialysis chamber as described previously.[47,48] Briefly, a 13 kDa MW cut-off dialysis membrane was used to separate the contents of the reaction vessels from a chamber with a continuously flowing buffer. 780 μL reaction mixtures were prepared similar to the mixtures described for the fluorometric enzyme-activity studies (4 fold the volumes) without the addition of a substrate. Activity control measurements were performed by the addition of pure DMSO instead of the inhibitor solutions. Due to the limited solubility of the inhibitors, the final inhibitor concentration was chosen to be 200 and 100 μM as 15% (v/v) DMSO was still well tolerated by Mpro (data not shown). The mixtures were transferred to the vessels of the dialysis chamber without incubation and dialyzed against a continuous flow of assay buffer containing 15% (v/v) DMSO (400 mL/h). Samples of 58.5 μL each were taken in triplicates at four different time points (0, 30, 60, 120 min), and reactions were initiated by the addition of 1.5 μL of Dabcyl-KTSAVLQSGFRKME-Edans (Genescript, New Jersey, USA) in a final concentration of 5 μM. Fluorescence was recorded over 10 min as described for fluorimetric enzyme activity studies.

Results and Discussion

Covalent Docking-Based VS Experiments

For the VS experiments, the X-ray structure of Mpro in complex with the potent covalent peptidomimetic inhibitor N3 (PDB 7BQY)[5] was used. This structure was chosen because it featured the highest resolution among the available ones when calculations were performed (1.70 Å). Analysis of the structure reveals that the binding site is located at the crevice between domains I and II of the enzyme where a Cys–His catalytic diad (C145 and H41) is present. N3 is lodged in this crevice in an extended conformation and forms a covalent bond through its Cβ vinyl group with C145 (Figure ). In this position, the lactam at position P1 (occupying the S1 site) forms a double H-bond with H163 and E166, while the N3 P2 Leu is inserted in a deep lipophilic cleft (S2). The P3 position is occupied by a solvent-exposed valine residue, while the alanine residue in P4 is engaging in van der Waals contacts with the protein at the S4 pocket. The terminal P5 and P1′ groups are partially exposed to the solvent although still able to form polar as well as van der Waals interactions with the S5 and S1′ pockets, respectively.

Figure 1

X-ray experimental complex between SARS-CoV-2 Mpro and the N3 inhibitor (PDB 7BQY).[5] The protein is represented as a green surface and sticks while the ligand as orange sticks. H-bond interactions are represented as dashed yellow lines.

Employing the above-described structure, covalent docking calculations were performed for all the ligands present in the in-house database. This chemical library can be roughly divided into two families of compounds. The first one (Figure ) is represented by a set of 37 dipeptidyl derivatives featuring a vinyl sulfone, a vinyl ester, two different vinyl ketones, and the vinyl nitrile as warheads. In all the derivatives, the homophenylalanine (hPhe) side chain is present at the P1 position, while aromatic and aliphatic amino acids occupy the P2 position. Finally, a set of aromatic and cyclic (namely, morpholine and piperazine) moieties are present at the P3 site.

Figure 2

Structural modifications of the dipeptidyl derivatives considered in VS calculations.

The second set of compounds which are non-peptidic ones, all feature the benzodiazepine core structure (60 compounds). These analogues were conceived as conformationally constrained peptidomimetics bearing a vinyl sulfone, a vinyl ester, a vinyl amide, a vinyl ketone, a vinyl phosphonate, and a vinyl nitrile as warheads (Figure ). The side chains of the hPhe, Gly, and Ile residues occupy the P1 position, while the P2 region is filled with variously substituted phenyl rings or a methyl group. High variability is present in the P3 position with aromatic, aliphatic, and positively charged substituents.

Figure 3

Structural modifications of the benzodiazepine derivatives considered in VS calculations.

Nine additional benzodiazepines were also available, featuring a vinyl ester, vinyl ketone, and vinyl nitrile warheads attached to the N1 atom of the benzodiazepine core through 1 to 3 carbon atoms (Figure ).

Figure 4

Benzodiazepine derivatives considered in VS calculations featuring different linkers between the core and the warhead.

These compounds were all docked into the selected three-dimensional (3D) structure of the Mpro enzyme, and results of these calculations were analyzed based on the ability of the docked ligands to recapitulate the N3/Mpro interaction pattern in the best ranking binding pose (i.e., the one having the lowest predicted binding free energy ΔGAD4). Such an analysis revealed that a number of the screened dipeptidyl derivatives are predicted to fit in the Mpro binding site sharing a docked pose that largely matches the co-crystal ligand (N3) interaction pattern (Figure ) with the P1, P2, and P3 groups being lodged in the S1, S2, and S3 pockets, respectively.

Figure 5

Schematic representation of the common structure of the dipeptidyl derivatives in the X-ray SARS-CoV-2 Mpro structure. The protein is represented as a green surface and sticks, while the modeled ligand is depicted as gray sticks. H-bond interactions are represented as dashed yellow lines. The warhead, P2, and P3 positions are represented as green, magenta, and red spheres, respectively. The P1 position (hPhe) is represented as blue sticks.

Regardless of the warhead, the peptides are predicted to form a covalent adduct with C145 and to form H-bonds with the residues lining the binding site cavity through their amide bonds (Figure ). Namely, the ligands P1 backbone NH contacts the backbone CO of H164, the P2 NH is in a favorable H-bonding position to interact with the side-chain CO of Q189, while the E166 backbone NH would form an H-bond with the P3 CO. The ligands P1 hPhe side chain is lodged in the enzyme S1 polar pocket, which comprises the residues F140, N142, E166, and H172. The P2 side chain inserts deeply into the hydrophobic S2 subsite, which consists of the side chains of H41, M49, and M165, as well as the alkyl portion of the side chain of D187. Compounds bearing a 4F-Phe P2 portion would project the partially negative fluorine substituent toward the S2 D187 side chain. It could be inferred that this repulsive electrostatic interaction might be detrimental to the binding process. Notably, a recent work by Zhang et al. reports that compounds with P2 moieties featuring a 4F-Phe are inactive against the SARS-CoV and HCoV-NL63 Mpro enzymes.[49] Indeed, the said proteases share a high degree of sequence identity with the Sars-CoV-2 Mpro. Thus, ligands featuring a 4F-Phe at the P2 site residue by engaging in unfavorable interactions with the D187 residue should not show interesting affinities against Mpro, even though this is not mirrored in the docking score, as the docking forcefield cannot adequately gauge this kind of interactions. Conversely, compounds that possess in P2 a methyl-cyclohexyl or a Leu side chain should provide a better fit for the S2 pocket. The P3 side chain is positioned in a pocket lined by the residues M165, L167, P168, Q189, and Q192. Given the pocket size, larger and bulkier substituents may be favored. As against rhodesain, several ligands bearing an electron-withdrawing group (EWG), mainly F atoms and NO2 groups on a phenyl or benzyl ring in P3, are supposed to enhance a π–π interaction with a nearby Phe. Indeed, the SARS-CoV-2 Mpro S3 pocket seems to include neither aromatic nor positively charged side chains that would establish charge–transfer interactions with an aromatic P3 moiety. However, it is not possible to rule out the possibility that EWGs could still be beneficial to the binding process, as they might enhance the π-stacking with surrounding π-faces of the residue backbones or coordinate water-mediated H-bonds with the enzyme. Finally, the majority of the ligands possess a vinyl ketone warhead group, as it was demonstrated to be the most efficient against rhodesain in our previous work.[26] In this latter study, through MD simulations, we computationally demonstrated that the increased ketone flexibility, if compared to the other tested warheads, was conducive to the formation of two H-bonds with the polar residues part of the S1′ site and contributed to the stabilization of the binding pose in the rhodesain active site. Interestingly, the attained docking solutions on Mpro should form an H-bonding interaction with the backbone NH of either G143 or C145 through a warhead acceptor atom. Thus, also in the Mpro case, the ketone flexibility might favor the formation of these polar contacts. On the other hand, the bulky benzenesulfonyl warhead might be well suited to occupy the rather wide S1′ pocket, which includes T25, T26, L27, and N142, while still retaining the capacity of engaging in an H-bond with the G143 or C145 backbone.

Docking results achieved for the 60 benzodiazepine derivatives were also analyzed to determine their propensity to form a covalent adduct featuring the interaction pattern established by the co-crystal N3 peptide in the Mpro enzyme. Unfortunately, none of the studied compounds was predicted to concurrently place the P1, P2, and P3 groups in the S1, S2, and S3 protein pockets. Thus, these ligands were not considered for further testing.

Compound Selection and Biological Evaluation

For some of the tested compounds, all belonging to dipeptidyl derivatives, AD4 was able to predict the lowest energy binding conformation in which the ligand can simultaneously make contacts with the S1, S2, and S3 pockets of the Mpro enzyme with its P1, P2, and P3 groups, respectively. In the present work, the existence of these specific ligand/protein interactions, rather than the ΔGAD4 values predicted through the AD4 native scoring function, was the sole criterion to select the most interesting candidates for subsequent biological evaluation. The Michael addition creates another chiral center in the β position to the warhead, and this leads to the formation of the (R)- and (S)-adducts. According to our calculations, the formation of both, the (R)- and (S)-complexes, should generally occur, as demonstrated by the similar ΔGAD4 obtained for the two complexes. Table reports the structures of the 15 compounds that were able to proficiently contact the enzyme (according to the selection criterion mentioned above) along with the predicted ΔGAD4 values of the two isomeric complexes. These compounds were selected for biological evaluation. Also, boceprevir and telaprevir (16 and 17, respectively) were used as positive controls being reported to inhibit the SARS-CoV-2 Mpro with IC50 values of 1.59 and 55.76 μM, respectively.[39]

For SARS-CoV Mpro, there is a wide range of KM values reported. This variability was attributed to the presence and position of a hexahistidine tag that could drastically impede the enzymatic activity.[50] Therefore, a construct with a native nsp5/nsp6 autocleavage site was designed, leading to a highly active, native Mpro. The KM of the nsp4/nsp5 cleavage site analogous substrate Dabcyl-KTSAVLQSGFRKME-Edans was determined to be 33.46 ± 3.01 μM using the GRAFIT software (Version 5.0.13; Erithacus Software Limited, East Grinstead, West Sussex, UK).

In an enzyme-inhibition assay, an initial screening of selected compounds was performed at a concentration of 100 μM and revealed that 6 and 10 inhibited SARS-CoV-2 Mpro more than 50% compared to the DMSO control (Figure A). For the two compounds, dilution series were prepared and reactions were initiated without prolonged incubation to determine their IC50 values to 47.2 ± 4.0 and 157.5 ± 9.3 μM, respectively.

Figure 6

Inhibitor screening and inhibition-mode determination in a fluorometric enzyme-inhibition assay. (A) Percentaged enzymatic activities of SARS-CoV-2 Mpro, treated with 100 μM inhibitor (compounds 1–15) or DMSO (ctrl. = 100%) in technical duplicates. (B) Relative activity of technical triplicates of samples from a dialysis assay at different time points. Error bars represent the respective standard errors.

To identify the inhibition mode of the most active compounds 6 and 10, a dialysis assay was performed, where a vessel, separated from a continuous flow cell by a semipermeable membrane, was filled with a reaction mix of the inhibitor and SARS-CoV-2 Mpro in the assay buffer. A reversibly inhibiting ligand would be expected to migrate through the membrane, and therefore, the enzymatic activity of the taken samples should recover over time (compared to a similarly treated DMSO control). In our experiment, however, no such increase in activity was observed over the monitored 120 min period. On the contrary, inhibition successively increased, indicating covalent reactions of the residual inhibitor with the enzyme (Figure B).

Binding-Mode Analysis of Compounds 6 and 10

Biological data demonstrated that compounds 6 and 10 inhibit the SARS-CoV-2 Mpro enzyme in the micromolar regimen. Interestingly, the two compounds feature:

The methyl vinyl ketone warhead.

An aliphatic amino acid in P2, namely, an l-cyclohexylalanine in compound 6 and an L-leu in 10.

A P3 para-substituted arene with an EWG, a fluorine atom in 6, and a NO2 group in 10.

These data support most of our initial inferences drawn from the docking experiments. It confirms the importance of hydrophobic interactions in the S2 pocket and the role played by the methyl vinyl ketone in the formation of a covalent adduct with C145 and possibly establishing an additional H-bond with G143 backbone NH (Figure ). Both ligands take favorable H-bonding contacts through their backbone atoms with residues lining the binding cavity (namely, G143, E166, and N189), while the P3 arene is lodged in the S3 pocket establishing van der Waals contacts, charge–transfer interactions with the π-faces of the nearby amino acids backbones, and H-bonds (compound 10 nitro group with N192 backbone NH).

Figure 7

Predicted binding mode of 6 (A) and 10 (B) into the X-ray SARS-CoV-2 Mpro structure. The protein is represented as a green surface and sticks while compounds 6 and 10 as pink and yellow sticks, respectively. H-bond interactions are represented as dashed yellow lines.

The above-described docking solutions would indicate a proper fit of these ligands into the Mpro enzyme and were instrumental for the experimental identification of two new inhibitors for this enzyme. Nevertheless, docking studies do not give information on the stability of the said interactions and the solvation effect. Thus, in the present work, MD simulations were attained on 6 and 10, the most active ligands, to probe the dynamic behavior of the ligand/protein interactions and give valuable hints for further structural modifications. To this end, the predicted 6/Mpro and 10/Mpro complexes were subjected to a 100 ns MD simulation to monitor the ligand/protein interactions along with the structural fluctuations through the sampled simulation time. The analysis was attained by examining the ligand root mean square deviations and fluctuations (L-rmsd and L-rmsf, respectively) to characterize the changes in the ligand atom positions.

As reported in Figure A,B, the docked poses calculated for 6 and 10 in the Mpro enzyme experience a different degree of stability. Most precisely, while 6 demonstrates to reach a very stable conformation (below 2 Å rmsd with respect to the starting frame), compound 10 experiences major fluctuations during the first 60 ns of the simulation time and reaches a more stable conformation during the last 40 ns of the production run.

Figure 8

rmsd (Å) plot over time (ns) of 6 (A) and 10 (B).

Figure A,B reports the main fluctuations for both 6 and 10 broken down by an atom. These plots were helpful in indicating how the different regions of the ligand interact with the protein and their entropic role in the recognition event. This analysis clearly indicates that the P1 position (hPhe) of both ligands is experiencing major movements while, differently from 6, the P3 position of 10 is also rather flexible. In general, comparison between 6 and 10 L-rmsf values obtained from this analysis demonstrated that 6 is able to establish more stable interactions with the enzyme if compared to 10.

Figure 9

L-rmsf plots broken down by an atom, corresponding to the reported two-dimensional structures of 6 and 10 (A,B, respectively).

Moreover, the nature of the ligand/protein interactions is different between 6 and 10. In particular, 6 establishes stable (>30% of the simulation time) H-bonds with the protein (Figure ). In particular, while confirming what already predicted by AD4 for the interactions with G143 and E166, MD simulations were also able to predict the formation of additional stable H-bonds with D155 and H164. Also, van der Waals contacts were found with I152 and P168 (9B and 9C). On the other hand, the afore described relocation of 10 (Figure B) during the MD simulation leads to a pose that is generally more solvent exposed than what predicted by AD4, in which water-bridged interactions with the protein take place as well as three H-bonds with H164, E166, and N189. All in all, the different stability of the ligand/protein interactions recorded for 6 and 10 would rationalize why the first is a more proficient Mpro inhibitor. Moreover, MD calculations allow us to infer that modifications of the hPhe P1 position on both ligands should be attempted to maximize the interactions with Mpro. In this respect, as depicted in Figure , N3 features a glutamine mimetic residue at the P1 position that engages two H-bonds with the side chains of residues lining the corresponding S1 pocket. Most probably, similar modifications of the most active compound 6 will result in higher inhibitory potencies against Mpro.

Figure 10

Protein interactions with compound 6 (A) and 10 (B) throughout the simulation categorized into: hydrogen bonds and hydrophobic, ionic, and water bridges. Panels C,D report a schematic of detailed 6 and 10 (C,D, respectively) atom interactions with the protein residues. Interactions that occur more than 30% of the simulation time in the selected trajectory (0.00 through 100.00 ns) are shown.

Conclusions

In the present work, covalent docking calculations were employed to virtually screen an in-house focused library of Michael acceptors in the pursuit of new SARS-CoV-2 Mpro lead compound inhibitors. Analysis of the docking results demonstrated that a number of the available dipeptidyl derivatives were predicted to form a covalent adduct with the reactive nucleophile thiol of C145 and concurrently occupy the S1, S2, and S3 pockets of the protein target. Enzyme inhibition of the most promising 15 compounds demonstrated that two of them (compounds 6 and 10) can inhibit the enzyme with IC50 values in the micromolar range. Besides, MD simulations achieved on the predicted 6/Mpro and 10/Mpro complexes allowed detecting the compound region (P1) that should first undergo future structural modifications to enhance its inhibition properties.

From the computational point of view, this paper demonstrates that, to identify active compounds against Mpro, the so-called “flexible side-chain method”, available within the AD4 software, can be successfully employed so that its application in the VS of larger ligand database of potential covalent binders can be envisaged. Moreover, the present account further underscores the custom nature of AD4 to address different docking issues in VS experiments such as the inclusion of protein flexibility,[51,52] the explicit consideration of the solvation effect,[53] the docking against nucleic acid targets,[54] as well as the formation of covalent adducts.[34] Even though our approach was ultimately successful, we recognize that a greater inclusion of protein flexibility in the screening might favor the identification of more hits potentially overseen in the rigid receptor docking campaign presented here. Certainly, an increasing number of Mpro crystal structures is now being released, and new campaigns could adopt an ensemble docking approach making use of parallel ranking.[52] As some of us have demonstrated, this protocol can lead to improved VS results provided that the considered protein structures feature a certain degree of binding site plasticity.[55] Also, we want to outline that in this work, the AD4 native scoring function was used to rank the different binding poses deriving from the docking of a single ligand and to pick the lowest energy binding conformation as a representative of all possible microscopic binding states. Rather, the presence of specific binding interactions in the docked compounds was the main criterion employed for the selections of the compounds to actually test. This was dictated by the fact that the well-known inaccuracies of the docking scoring functions are further burdened, in covalent docking, by the neglecting of the whole contribution deriving from the covalent bond formation.