key: cord-0649224-9m47cbwf
authors: Menendez, Cintia A.; Bylehn, Fabian; Perez-Lemus, Gustavo R.; Alvarado, Walter; Pablo, Juan J. de
title: Molecular Characterization of Ebselen Binding Activity to SARS-CoV-2 Main Protease
date: 2020-05-20
journal: nan
DOI: nan
sha: efd0db816c414f8900dac9b53bc6a89da90fd00a
doc_id: 649224
cord_uid: 9m47cbwf

The Coronavirus Disease (COVID-19) pandemic caused by the SARS-coronavirus 2 (SARS-CoV-2) urgently calls for the design of drugs directed against this new virus. Given its essential role in proteolytic processing, the main protease Mpro has been identified as an attractive candidate for drugs against SARS-CoV-2 and similar coronaviruses. Recent high-throughput screening studies have identified a set of existing, small-molecule drugs as potent Mpro inhibitors. Amongst these, Ebselen (2-Phenyl-1,2-benzoselenazol-3-one), a glutathione peroxidase mimetic seleno-organic compound, is particularly attractive. Recent experiments suggest that its effectiveness is higher than that of other molecules that also act at the enzyme catalytic site. By relying on extensive simulations with all-atom models, in this study we examine at a molecular level the potential of Ebselen to decrease Mpro catalytic activity. Our results indicate that Ebselen exhibits a distinct affinity for the catalytic site cavity of Mpro. In addition, our molecular models reveal a second, previously unkown binding site for Ebselen in the dimerization region localized between the II and III domains of the protein. A detailed analysis of the free energy of binding indicates that the affinity of Ebselen to this second binding site is in fact significantly larger than that to the catalytic site. A strain analysis indicates that Ebselen bound between the II-III domains exerts a pronounced allosteric effect that regulates catalytic site access through surface loop interactions, and induces a displacement and reconfiguration of water hotspots, including the catalytic water, which could interfere with normal enzymatic function. Taken together, these findings provide a framework for the future design of more potent and specific Mpro inhibitors, based on the Ebselen scaffold, that could lead to new therapeutic strategies for COVID-19.

, hydrogen bonding interactions are formed between the carbonyl oxygen of Ebselen and the Asn142 and Gln189 side chains. One can also appreciate the hydrophobic contacts between Ebselen and Met165, Pro168, Met49, His164 and both residues belonging to the catalytic dyad, His41 and Cys145. Figure 2B shows a representative configuration of Ebselen at the intersection between domains II and III. There are distinct hydrophobic contacts with Phe294, Pro108, Ile200, Val202, His246, Thr292, Ile249, Pro132 and Ile249 (some residues are not shown for clarity). Highly dynamic hydrogen bonds were observed between Ebselen and the Gln107, Gln110, and Hie246 side chains. The results shown in Figure 1 were generated on the basis of direct molecular dynamics simulations. In order to arrive at a quantitative estimate of the binding affinity of each site, we used thermodynamic integration to determine the absolute binding free energy for the catalytic site and the Domain II-II site. More specifically, we used Particle Mesh Ewald Molecular Dynamics (PMEMD) as implemented in Amber 18 (33) with 11 windows per integration and 10ns per window. In addition, multiple runs starting from the most probably binding cluster identified previous MD simulations were considered. In this way, three independent replicas for each site were used to calculate averages. The results are shown in Table 1 Consistent with the probability density maps shown in Figure 1 , the absolute binding free energies corresponding to both sites are negative, serving to underscore the thermodynamic stability of the Ebselen-Mpro complex at both sites. Here we note the binding affinity of Ebselen at the second binding site, located between the II and III domains, is in fact greater that at the catalytic site. A previous report suggested that this small molecule could also inhibit Mpro through non-covalent binding, particularly because it has been found to exhibit a stronger inhibition effect than other compounds that were also able to modify Cys145 in the catalytic dyad quantitatively (5) . The results presented here could explain Ebselen's high enzymatic inhibition, even when this compound could only partially modify Cysteine 145.

Having identified a distant binding site for Ebselen between domains II and III, we examined the role of this binding site, if any, on the catalytic site of the protein.

To do so, we determined the local strain induced by the binding throughout the protein.

To examine potential synergies of drug binding to both sites, we also considered a third scenario in which both sites (catalytic and distant) are simultaneously occupied by Ebselen. Note that local strain provides a measure of local deformation that filters out non-trivial, functional conformational changes from non-functional ones, and it is therefore ideally suited to highlight how a perturbation (i.e. Ebselen binding) at one site induces conformational changes at other, potentially distant, sites. Note that different measures of strain have recently been used in studies of allostery in proteins (17, 18, 19) . Here, we used a method originally introduced by our group in the context of local strain in polymeric glasses. From the β-factor analysis, it is evident that when Ebselen is bound to the Domain-II-III interface ( Figure 3A : blue lines), the dynamics of the 44-52 loop, which flank the catalytic site, is significantly altered. On the other hand, when Ebselen is bound at both sites simultaneously there is a clear global reduction of the thermal fluctuations among the receptor except for one single region located around residues 137-140, where enhanced flexibility is apparent. Apart from these residues, Ebselen bound at both sites shows the lowest β-factor values, however, at this specific point, this system exhibits as high flexibility as Mpro-apo protein ( Figure 3A : green line). The strain analysis helps disentangle the effects of functional and non-functional fluctuations and provides a more detailed view of the effects of Ebselen binding. As seen in Figure 3A for the β-factor analysis, when Ebselen is bound to the domain interface it produces a large strain at the 44-52 loop ( Figure 3B : blue line). Likewise, the 185-201 loop that also flanks the catalytic site, as well as the region comprising residues 137-140, exhibit a high strain, which is not immediately apparent in the β-factor. When both molecules are bound simultaneously to Mpro ( Figure 3B , black line), a high strain signal is also exhibited around residues 137-140, as expected based on the aforementioned β-factor results. Finally, when Ebselen is located in the catalytic site, a much lower strain is showed around all these aforementioned three regions ( Figure 3B , red line). Figure 3C shows that the strain is primarily localized at the two loops flanking the catalytic site (44-52 and 185-201 loops), as well as a loop in catalytic site (residues 22-25). A highly strained region also appears around residues 137-140. To understand why the 137-140 residues show such a large strain when Ebselen is bound to both sites we turn our attention to the molecular images shown in Figure 4 . In the close-up of the residues shown in Figure 4A , a specific backbone hydrogen bond is formed between Lys137 and Phe140 due to the binding of Ebselen between domains II and III; this conformational change induced by the presence of Ebselen causes the high strain shown in Figure 3B . This conformational change takes place in the middle region between the catalytic site (Cys145) and the binding cleft between domains II and III (Pro132), which points to the relevance of these residues.

From these results, it is evident that when Ebselen binds between the II-III domains, it exerts a pronounced allosteric effect that affects the loops that regulate access to the catalytic site (44-52 and 185-201 loops). In addition, it affects the residues 137-140, where a specific backbone hydrogen bond is formed between Lys137 and Phe140. The exact role of this conformational change will be the subject of future studies, but the results presented here show that it acts as a relay between domain III and the catalytic site.

Given that conformational changes in the catalytic site are observed when Ebselen binds to Mpro far away from this specific region, we turn our attention to the hydration characteristics of the catalytic site in the Mpro-Ebselen complex (when bound to domain II-III interface) and in the Mpro-apo structure (PDB code: 6m03) (20) . AQUADUCT 1.0.5 (21) was used to analyze the water structure and water flux in the Mpro protein, with a time window of 50 ns, and sampling every 1 ps. The results are shown in Figure 5 . ) shows that Ebselen binding between Domain II-III leads to fewer water inlets compared with the Mpro apo state; this is indicative of a water flux reduction in the catalytic site due upon Ebselen binding, even though it happens far from the catalytic site. From the lower panel in Figure 5 , it is evident that a volume reduction of the catalytic site occurs when Ebselen is located between domain II-III compared with the apo state; the Maximum Available Volume (MAV) for water in the catalytic region is around 50% smaller in this case. Note that in the apo protein there is a catalytic water that forms a catalytic triad together with Cys145 and His41 (16, 23) , and this catalytic water (red in Figure 5A lower panel) is preserved and remains close to His41 in the simulations. In contrast, it is clear from Figure 5B (lower panel) that the presence of Ebselen induces a displacement and reconfiguration of water hotspots, including the catalytic water (red). Importantly, these effects could prevent the normal enzymatic function of Mpro, as this catalytic water displacement might damage the catalytic triad that is required for protein activity (24) . Similarly, the observed pocket volume reduction might affect the accessibility of the polyprotein that Mpro cleaves, thereby reducing enzymatic function.

As mentioned before, Mpro is an attractive drug target against the COVID-19 virus due to is its central role in the viral life cycle. A previous structural and evolutionary investigation suggested that the SARS-CoV-2 Mpro is not a suitable target for de novo development of inhibitors or the repurposing of drugs against the previous SARS coronavirus (16) . That study, however, only compared the active sites in Mpro for COVID-19 and the highly similar previous SARS-CoV Mpro in terms of flexibility and plasticity, where major differences in both shape and size were observed indicating that repurposing SARS drugs for COVID-19 may not be effective. Based on their evolutionary analysis, these authors also pointed out that the virus's mutability will pose further challenges to treatments against the COVID-19 Mpro protein. An alternative to this discouraging scenario, however, would be to target the region between the II and III domains, which is implicated in dimer formation.

Here we find that there are two, highly probable interaction sites between SARS-CoV-2 Mpro and Ebselen. One is located within the catalytic cavity and, importantly, a second site is in the region between the II and III domains, which is essential for Mpro dimerization (15) . Detailed calculations of the free energy reveal a higher binding affinity of Ebselen to the Domain II-III than to the catalytic site.

The strain analysis reveals that Ebselen bound between the II-III domains exert a pronounced allosteric effect that affects the loops regulating access to the catalytic site. In addition, it also affects residues 137-140, where a specific backbone hydrogen bond is formed between Lys137 and Phe140.

The catalytic site water analysis indicates that the proposed allosteric inhibition by Ebselen could occur through a volume reduction of the catalytic pocket, and a reconfiguration of water hotspots in that region. Given the catalytic role of water in this enzyme's activity, these effects could act to prevent the regular enzymatic function of the SARS-CoV-2 Mpro protein. A previous evolutionary study (16) showed that the mutation of a few residues belonging to the catalytic site is energetically unfavorable. Therefore, residues such as P39, R40, P52, G143, G146, or L167 could be considered as key anchoring residues for Mpro inhibitor design. The insights put forth in this work have the potential to facilitate the rational molecular design of new analogs, based on the Ebselen scaffold, that result in anchoring at those positions. This alternative approach, uncovered through the extensive molecular simulations presented here, highlights the need to develop reliable, high-throughput methods to screen drug-protein interactions at the molecular level and that incorporate the role of explicit water morelcules. We conclude by pointing out that our focus in this work has been on non-covalent complexes between Ebselen-Mpro. In a subsequent stage, we also plan to investigate covalent complexes involving Cys145.

A total of more than 6 μs of classical MD simulations of SARS-CoV-2 Mpro apo state and Mpro-Ebselen complex were run using AMBER18 (28) simulation package (3 μs: ebselen as a molecular probe; 2.4 μs: shear strain analysis; 100ns: water structure and flux analysis; 990ns: free energy analysis). The receptor initial configuration for the Mpro-Ebselen system was taken from recently reported structure for Mpro-N3 inhibitor (5) (PDB ID: 6lu7); where the inhibitor and crystallographic water molecules were removed before starting simulations. Force field parameters for Ebselen were determined using the Antechamber program and described by the General Amber Force Field (GAFF) (26, 27) . The partial atomic charges were determined by the restrained electrostatic potential (RESP) fitting technique. Those electrostatic potential calculations were performed at the HF/6-31G level with Gaussian 09. For Mpro apo simulations we used the recently reported structure (20) , (PDB ID: 6m03). Approximately 20,000 TIP3P water model molecules and 4 Na + ions were added. All simulations were carried out using the ff14SB force field (29) . The simulation protocol included a first minimization of 7000 steps, involving 3500 steepest descent steps followed by 3500 steps of conjugate gradient energy minimization, where constraints were applied on the protein heavy atoms (force constant 500 kcal x mol -1 x Å 2 ) and a second minimization (7000 steps) with no constraints of conjugate gradient energy minimization. Next, during the first equilibration, the temperature was gradually increased from 0 K to 300 K over 50 ps using a Langevin thermostat with a temperature coupling constant of 1.0 ps in the canonical ensemble. Density equilibration and production runs were carried out using a constant pressure ensemble (NPT). All simulations were performed using periodic boundary conditions and a 2 fs time step. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald method with a non-bonded cut-off of 10 Å and the SHAKE algorithm was used to implement rigid constraints.

In order to determine the most probable interaction sites between Mpro-Ebselen (global hot-spot), 15 replicas of 200 ns each were run (a total of 3 μs), where configurations were saved every 20 ps. The CPPTRAJ (30) software was used to process the trajectories, where at the first stage all trajectories were align by means of minimization the distance among protein backbone atoms (C, N, CA). Grid command was used to track Ebselen molecule and to produce number density map, where the grid resolution was selected to be 0.5 Å.

Initial Mpro-Ebselen complex configurations were selected from previous global binding affinity study, with Ebselen being bound to the catalytic site, to the intersection between Domain II-III, and to both of these sites simultaneously. Three replicas of 200 ns for each of these initial setups were run. As a reference simulation, we used the Mpro-apo structure (PDB code: 6m03), and 3 replicas of 200 ns each were run. In this way, a total of 2400 ns of classical MD were run.

To apply the strain formalism from continuum theory to discrete, atomistic systems, differential operators replace the derivatives (17, 18, 19, 31, 32) . A radius around each central atom containing other atoms defines the local neighborhood around the central atom. The instantaneous position of atom at any timestep in the MD is , and the position of the same atom at any timestep of the reference simulation is 0, . To first order, the distances between atom and its neighbors are related through the deformation matrix F by − = ( 0, − 0, ), which forms an overdetermined system of linear equations, and an optimized * is sought by minimizing the difference between the actual distances and the projected distances to an affine deformation:

The atomic strain tensor is then found by = since proteins are generally incompressible (17, 19) . For this analysis, a radius of 10 Å around each atom is considered (19) , and only the Cα atoms are used in the calculations. The reference simulation is the apo protein trajectory and the strain is then measured using the Ebselen simulations to elucidate the effect of binding of Ebselen at the different sites.

In addition, β-factors were estimated over the same trajectories using the atomicfluct command of the CPPTRAJ module of Amber18.

To analyze the water structure and water flux in the Mpro protein of SARS-CoV-2 we used AQUADUCT 1.0.5 (21) . We obtained Inlets, Maximum Available Volume (MAV) and Hotspots for water molecules in the Mpro catalytic region. For calculations, this region was defined as a 5 Å sphere around the center of geometry of the active site residues (H41, C145, H164, D187) (16) . The Mpro protein was studied in two different scenarios, with Ebselen in the Domain II-III site and the apo protein with no ligand. The time window used in both calculations was 50ns, sampling every 1 ps. Images were created with open-source PyMOL (22) .

The absolute binding free energy is defined as: ΔGbinding=ΔG L -ΔG RL , where ΔG RL is the free energy change of Ebselen annihilation in the Mpro complex, and ΔG L is the free energy change of Ebselen annihilation in water. To calculate these free energy changes, we use Thermodynamic Integration (TI) implemented in PMEMD for Amber 18. We use the one step anhihilation protocol with soft core potentials (34) . In addition, we adopted a simple approach with multiple runs starting from the most probably binding cluster estimated from previous MD simulations. In this way, three independent replicas for each site were taken into account, as well as three replicas for Ebselen solvated in pure water. Eleven equally spaced windows were used (ΔLambda=0.1) with 10ns of simulation time per window. To keep the ligand from wandering in TI calculations, we used a soft restraint of 10 kcal/molÅ 2 (35) .

A novel coronavirus from patients with pneumonia in China

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

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

A new coronavirus associated with human respiratory disease in China

Structure of M pro from COVID-19 virus and discovery of its inhibitors. bioRxiv preprint

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

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

A safe lithium mimetic for bipolar disorder

Development of ebselen, a glutathione peroxidase mimic, for the prevention and treatment of noiseinduced hearing loss

Safety and efficacy of ebselen for the prevention of noise-induced hearing loss: a randomised, double-blind, placebo-controlled, phase 2 trial. The Lancet

Repurposing ebselen for treatment of multidrug-resistant staphylococcal infections

Synergistic antibacterial effect of silver and ebselen against multidrug-resistant Gram-negative bacterial infections

2-Phenyl-1, 2-benzisoselenazol-3 (2H)-one containing pharmaceutical preparations and process for the treatment of rheumatic diseases

Effects of the potential lithium-mimetic, ebselen, on impulsivity and emotional processing

Maturation Mechanism of Severe Acute Respiratory Syndrome (SARS) Coronavirus 3C-like Proteinase

Molecular Dynamics Simulations Indicate the COVID-19 Mpro Is Not a Viable Target for Small-Molecule Inhibitors Design

Strain analysis of protein structures and low dimensionality of mechanical allosteric couplings

A deformation gradient tensor and strain tensors for atomistic simulations

Colloquium: Proteins: The physics of amorphous evolving matter

The crystal structure of COVID-19 main protease in apo form

AQUA-DUCT 1.0: structural and functional analysis of macromolecules from an intramolecular voids perspective

Pymol: An open-source molecular graphics tool. CCP4 Newsletter On Protein Crystallography

Coronavirus Main Proteinase (3CLpro) Structure: Basis for Design of Anti-SARS Drugs

The role of conserved water molecules in the catalytic domain of protein kinases

Virtual screening and repurposing of FDA approved drugs against COVID-19 main protease

Automatic atom type and bond type perception in molecular mechanical calculations

Development and testing of a general AMBER force field

ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data

On identifying collective displacements in apo-proteins that reveal eventual binding pathways

Dynamics of viscoplastic deformation in amorphous solids

Soft-core potentials in thermodynamic integration: Comparing one-and twostep transformations

Absolute Binding Free Energy Calculation and Design of a Subnanomolar Inhibitor of Phosphodiesterase-10