key: cord-312664-tgpaidhp
authors: Liang, Julia; Pitsillou, Eleni; Karagiannis, Chris; Darmawan, Kevion K; Ng, Ken; Hung, Andrew; Karagiannis, Tom C.
title: Interaction of the prototypical α-ketoamide inhibitor with the SARS-CoV-2 main protease active site in silico: Molecular dynamic simulations highlight the stability of the ligand-protein complex
date: 2020-05-28
journal: Comput Biol Chem
DOI: 10.1016/j.compbiolchem.2020.107292
sha: 
doc_id: 312664
cord_uid: tgpaidhp

The severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) causes an illness known as COVID-19, which has been declared a global pandemic with over 2 million confirmed cases and 137,000 deaths in 185 countries and regions at the time of writing (16 April 2020), over a quarter of these cases being in the United States. In the absence of a vaccine, or an approved effective therapeutic, there is an intense interest in repositioning available drugs or designing small molecule antivirals. In this context, in silico modelling has proven to be an invaluable tool. An important target is the SARS-CoV-2 main protease (M(pro)), involved in processing translated viral proteins. Peptidomimetic α-ketoamides represent prototypical inhibitors of M(pro). A recent attempt at designing a compound with enhanced pharmacokinetic properties has resulted in the synthesis and evaluation of the α-ketoamide 13b analogue. Here, we performed molecular docking and molecular dynamics simulations to further characterize the interaction of α-ketoamide 13b with the active site of the SARS-CoV-2 M(pro). We included the widely used antibiotic, amoxicillin, for comparison. Our findings indicate that α-ketoamide 13b binds more tightly (predicted GlideScore = -8.7 and -9.2 kcal/mol for protomers A and B, respectively), to the protease active site compared to amoxicillin (-5.0 and -4.8 kcal/mol). Further, molecular dynamics simulations highlight the stability of the interaction of the α-ketoamide 13b ligand with the SARS-CoV-2 M(pro) (ΔG = -25.2 and -22.3 kcal/mol for protomers A and B). In contrast, amoxicillin interacts unfavourably with the protease (ΔG = +32.8 kcal/mol for protomer A), with unbinding events observed in several independent simulations. Overall, our findings are consistent with those previously observed, and highlight the need to further explore the α-ketoamides as potential antivirals for this ongoing COVID-19 pandemic.

effective therapeutic, there is an intense interest in repositioning available drugs or designing small molecule antivirals. In this context, in silico modelling has proven to be an invaluable tool. An important target is the SARS-CoV-2 main protease (M pro ), involved in processing translated viral proteins. Peptidomimetic α-ketoamides represent prototypical inhibitors of M pro . A recent attempt at designing a compound with enhanced pharmacokinetic properties has resulted in the synthesis and evaluation of the α-ketoamide 13b analogue. Here, we performed molecular docking and molecular dynamics simulations to further characterize the interaction of α-ketoamide 13b with the active site of the SARS-CoV-2 M pro . We included the widely used antibiotic, amoxicillin, for comparison. Our findings indicate that α-ketoamide 13b binds more tightly (predicted GlideScore = -8.7 and -9.2 kcal/mol for protomers A and B, respectively), to the protease active site compared to amoxicillin (-5.0 and -4.8 kcal/mol). (1) . A novel coronavirus was isolated and designated severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), causing coronavirus disease 2019 . As of April 16, 2020 , this ongoing global health emergency has resulted in over 2,000,000 confirmed cases in 185 countries and regions, with more than 25% of confirmed cases in the United States (2) . The global mortality rate has been estimated to be 5.7%, with higher mortality occurring among the elderly (3) . The majority of deaths have occurred among adults aged greater than 60 years and those with serious underlying health conditions, with the highest fatality in those aged greater than 85 years ranging from 10% to 27% in the United States (4, 5) . Differences in disease prevalence are affected by sex, with data indicating that there is a higher prevalence of COVID-19 among men (6, 7) . The majority of early cases were linked to exposure to the Huanan Seafood Wholesale Market, potentially through zoonotic transmission (8) . Human-to-human transmission of SARS-CoV-2 was subsequently found to occur, with an attack rate within families of 83% suggestive of its high transmissibility (9, 10).

The current outbreak of SARS-CoV-2 follows that of recent outbreaks of severe acute respiratory syndrome coronavirus (SARS-CoV) in 2002 and the Middle East respiratory syndrome coronavirus (MERS-CoV) in 2012 (11) . These coronaviruses are both zoonotic pathogens, with bats serving as the primary reservoir (12) . Masked palm civets were the intermediate reservoir for SARS-CoV, and dromedary camels for MERS-COV, where zoonotic transmission to humans subsequently occurred (12) . While SARS-CoV-2 appears to have lower fatality rates than SARS-CoV (9.5%) and MERS-CoV (34.4%), it has a greater ability to spread (11, 13) . Like SARS-CoV, the pathogenesis of SARS-CoV-2 involves the binding of its spike protein to angiotensin converting enzyme-2 (ACE2) in the host (14, 15) . When cleavage occurs between the S1 and S2 subunits, the spike protein becomes activated for J o u r n a l P r e -p r o o f membrane fusion for entry into the host cell (14, 15) . ACE2 is expressed on numerous tissues in the nasopharynx and intestinal epithelia, particularly in type II alveolar cells in the lung (16) (17) (18) . Following entry of the virus into the host cells, viral RNA attaches to the host ribosome for translation of large polyproteins that are processed via proteolysis into components for new virions (19, 20) . Along with the papain-like protease, the coronavirus main protease (M pro ) is responsible for this proteolysis (19) . Encoded by open reading frame 1 (ORF1) of the genome as non-structural protein 5 (Nsp5), M pro cleaves at 11 sites in the polyproteins (19) .

To date, there is an absence of a vaccine and a lack of effective antiviral therapeutics against SARS-CoV-2. Therefore, there is an intense interest in identifying compounds that may interact with key viral molecular targets. Due to their functional importance and high degree of conservation among coronaviruses, M pro s have become an important target in the design of anti-coronaviral drugs (19, 21) . The structure of the SARS-CoV-2 M pro was initially solved by Jin et al. in late January of this year (22) , accelerating the search for drugs that may act as lead compounds. Following the 2002 SARS outbreak, work by Hilgenfeld at al. aimed at designing compounds with broad-spectrum anti-coronaviral activity, focussing on main proteases (19, 23) . Previously, they found that peptidomimetic α-ketoamides were potential candidates for broad-spectrum inhibitors of coronavirus and enterovirus replication (24) . Most recently, work aimed at improving the biological properties to produce an inhibitor specific for the SARS-CoV-2 M pro resulted in the potential antiviral agent α-ketoamide 13b (25) . It was found that compound 13b demonstrates binds to the substrate-binding cleft and exhibits antiviral activity in vitro, inhibiting the SARS-CoV-2 M pro with IC50 = 0.67 ± 0.18 µM (25) .

Here, our aim was to further investigate the interaction of the α-ketoamide 13b with the SARS-CoV-2 M pro in silico. To highlight the importance of molecular dynamics simulations, we J o u r n a l P r e -p r o o f compared the properties of α-ketoamide 13b, with one of the most widely prescribed antibiotics, amoxicillin. Amoxicillin was chosen for comparison for two reasons: 1) although it does not possess antiviral properties, it remains a mainstay as a frontline therapy for viral infections, including COVID-19, presumably to protect from opportunistic secondary bacterial infections, and 2) our initial screening indicates that it binds to the active site of the SARS-CoV-2 M pro akin to -ketoamide 13b, albeit with lower affinity.

Preparation of systems and docking calculations were carried out using the Schrodinger Suite Crystallographic water molecules were removed. The homodimer protein complex was prepared using the Protein Preparation Wizard (28) . This was used to assign bond orders, add hydrogens, create zero-order bonds to metals, and create disulphide bonds. Hydrogen bonds were assigned and optimised, followed by restrained energy minimization. Ligand structures were pre-processed using LigPrep (29) , for the generation of ionization and stereoisomer variants of input molecules to obtain structures with optimised geometry.

Two receptor grids of 20 x 20 x 20 Å in size were generated around the active site of the protease, centroid to residues GLY-23, THR-24, GLY-143, HIS-163, THR-190, and ALA-191 on each chain. Ligands were docked to each chain separately. Docking was carried out using the Quantum Mechanics-Polarized Ligand Docking (QPLD) workflow (30) of Schrodinger.

Initial docking was performed using the extra precision (XP) scoring function of Glide (31) .

Partial charges on ligand atoms were then calculated using quantum mechanical methods using J o u r n a l P r e -p r o o f the 'accurate' setting in Jaguar (32) . Ligands were re-docked using the calculated charges with XP docking mode of Glide, and the final pose was selected based on GlideScore.

Classical MD simulations were performed using GROMACS 2018.2 software(33, 34) with the CHARMM27 force field (35, 36) . Ligand topology was generated using SwissParam (37) .

Protein-ligand complexes were solved using TIP3P water (38) in a dodecahedral box with a minimum of 2.0 nm distance between any protein atom to the closest box edge. Sodium ions were added to the solvated system to neutralise the charge. Energy minimisation was performed using a steepest-descent gradient method for a maximum of 50,000 steps. Each complex was then restrained using an isothermal-isochloric (NVT) ensemble and isothermal-isobaric ensemble (NPT) for 100 ps. Temperature was maintained at 310 K with a modified Berendsen thermostat (39) , and pressure at 1.0 bar with the Parrinello-Rahman barostat (40) . Bond lengths were constrained using the LINCS algorithm (41) , with long-range electrostatic forces calculated using the particle-mesh Ewald scheme (PME) (42) (grid spacing 0.16 nm). Cutoff ratios of 1.2 nm for Coulomb and van der Waals potentials were used for the calculation of short-range nonbonded interactions. Simulations were carried out for 100 ns with a time-step of 2 fs in triplicate, with random generation of velocities according to a Maxwell distribution.

Visual Molecular Dynamics 1.9.3 (43) was utilised for analysis and visualisation of trajectories.

Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) was utilised for the quantification of free energy calculations (44) . This was performed using the g_mmpbsa tool (45) . MM-PBSA calculations were performed on 1 ns segments of the triplicate stabilised trajectories (46) . Energy contributions from electrostatic, van der Waals, and polar solvation terms were calculated using the adaptive Poisson-Boltzmann Solver (APBS) (47) . Grid spacing was set to 0.05 nm, and values of 80 and 2 were used for solvent dielectric constant and solute J o u r n a l P r e -p r o o f dielectric constant, respectively. Solvent-accessible surface area (SASA) was used to approximate the non-polar energy contribution, with the probe radius set to 0.14 nm. Entropic energy terms were excluded from the calculations.

The α-ketoamide 13b ligand binds with relatively high affinity to the active site of the SARS-

Docking was performed using the QPLD workflow of Schrodinger to obtain a rigorous estimate of ligand binding affinities to the active site. The receptor grid was centred around the active site residues for a comprehensive search of binding poses using partial atomic charges calculated using quantum mechanical methods. The α-ketoamide 13b ligand bound to the protease with a GlideScore of -8.8 kcal/mol to the monomer, compared to -5.2 kcal/mol for amoxicillin. For binding to the dimer structure of the protease (Table 2) 

It is noted that from visual inspection, while α-ketoamide 13b remains strongly bound to the active site of the protease throughout the trajectory (Media S1), amoxicillin does not remain Table 2 indicates that van der Waals interactions are the predominant driving force for binding of α- binding. Polar solvation energy was also significantly more unfavourable compared to binding with α-ketoamide 13b. While MM-PBSA is a commonly used method for estimating free energy, it should be noted that entropy terms are excluded from calculations. While MM-PBSA methods have been used to rescore poses obtained through molecular docking, it should be noted that values produced should be treated as relative differences across a set of ligands, rather than absolute (48, 49). Nevertheless, the results presented here demonstrate that αketoamide 13b is able to bind strongly to the substrate binding site of the protease.

Residue energy contributions are decomposed in Figure 4 , where the average energy contributions are shown on a per-residue basis for the entire protease. For α-ketoamide 13b, energy contributions are largely confined to the protomer on which they are bound. Residue energy contributions are mostly favourable, shown as negative peaks below the x-axis.

The residue with the most negative energy contribution corresponding to the most favourable interaction is MET-49 on both protomers. This residue lies within the region where the large fluctuations in RMSF were observed in Figure 2 . A more favourable energy contribution was observed for this residue on protomer A on the dimer of -2.40 kcal/mol compared to -1.26 kcal/mol for protomer B. This is in line with the slightly stronger ΔG of for α-ketoamide 13b binding to protomer A compared to protomer B (Table 2) . MET-49 is located in the S2 subsite of the substrate binding pocket, previously found to form a 'lid' in the closely related SARS-CoV M pro S2 site facilitating hydrophobic interactions that enable inhibition of the enzyme with α-ketoamides (24) . Another residue contributing favourably to binding to both protomers where they found that GLU-166 adopted an inactive conformation in protomer B (25) . GLU-166 is a key residue essential for catalytic activity, interacting with NH2 terminal residues of each protomer at the dimer interface, shaping the S1 pocket of the substrate binding site (25, 52) . 

TCK, KN, and AH conceptualized the aims and methodology and were involved in supervision. JL performed formal data analysis, was involved in data curation, and produced the first draft of the manuscript. CK performed formal data analysis and was involved in data curation. EP and KKD performed formal data analysis and validation. All authors contributed to editing and reviewing the manuscript. Amoxicillin -5.0 -4.8 

Situation Report

An interactive web-based dashboard to track COVID-19 in real time

Real estimates of mortality following COVID-19 infection. The Lancet Infectious Diseases

Severe outcomes among patients with coronavirus disease 2019 (COVID-19)-United States

Novel Coronavirus Pneumonia Emergency Response Epidemiology Team. The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases

Zhonghua liu xing bing xue za zhi= Zhonghua liuxingbingxue zazhi

Sex difference and smoking predisposition in patients with COVID-19. The Lancet Respiratory Medicine

A novel coronavirus outbreak of global health concern. The Lancet

Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus-Infected Pneumonia

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster

SARS-CoV-2 and COVID-19: The most important research questions

A Novel Coronavirus Emerging in China -Key Questions for Impact Assessment. The New England journal of medicine

SARS and MERS: recent insights into emerging coronaviruses

The many estimates of the COVID-19 case fatality rate. The Lancet Infectious Diseases

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

Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein

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

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

High expression of ACE2 receptor of 2019-nCoV on the epithelial cells of oral mucosa

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

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

Structures of two coronavirus main proteases: implications for substrate binding and antiviral drug design

Structure of Mpro from COVID-19 virus and discovery of its inhibitors

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

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 alpha-ketoamide inhibitors

Protein interfaces, surfaces and assemblies service PISA at European Bioinformatics Institute

Identifying and characterizing binding sites and assessing druggability

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

Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach

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

Berendsen HJC, van der Spoel D, van Drunen R. GROMACS: A message-passing parallel molecular dynamics implementation

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

Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models

A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields

SwissParam: a fast force field generation tool for small organic molecules

Comparison of simple potential functions for simulating liquid water

Molecular dynamics with coupling to an external bath

Crystal Structure and Pair Potentials: A Molecular-Dynamics Study

LINCS: a linear constraint solver for molecular simulations

Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems

VMD: visual molecular dynamics

Electrostatics of nanosystems: Application to microtubules and the ribosome

g_mmpbsa-A GROMACS Tool for High-Throughput MM-PBSA Calculations

Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations

Brown SP, Muchmore SW. Large-Scale Application of High-Throughput Molecular Mechanics with Poisson−Boltzmann Surface Area for Routine Physics-Based Scoring of Protein−Ligand Complexes

Investigation of MM-PBSA rescoring of docking poses

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

The papain-like protease of severe acute respiratory syndrome coronavirus has deubiquitinating activity

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

We