key: cord-0795703-emzr3dzb
authors: Menéndez, Cintia A.; Byléhn, Fabian; Perez-Lemus, Gustavo R.; Alvarado, Walter; de Pablo, Juan J.
title: Molecular characterization of ebselen binding activity to SARS-CoV-2 main protease
date: 2020-09-11
journal: Sci Adv
DOI: 10.1126/sciadv.abd0345
sha: 0f4d23e9e9e98d1c01dfe023066963db75d97ece
doc_id: 795703
cord_uid: emzr3dzb

There is an urgent need to repurpose drugs against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Recent computational-experimental screenings have identified several existing drugs that could serve as effective inhibitors of the virus’ main protease, M(pro), which is involved in gene expression and replication. Among these, ebselen (2-phenyl-1,2-benzoselenazol-3-one) appears to be particularly promising. Here, we examine, at a molecular level, the potential of ebselen to decrease M(pro) activity. We find that it exhibits a distinct affinity for the catalytic region. Our results reveal a higher-affinity, previously unknown binding site localized between the II and III domains of the protein. A detailed strain analysis indicates that, on such a site, ebselen exerts a pronounced allosteric effect that regulates catalytic site access through surface-loop interactions, thereby inducing a reconfiguration of water hotspots. Together, these findings highlight the promise of ebselen as a repurposed drug against SARS-CoV-2.

A coronavirus of zoonotic origin, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is the etiological agent responsible for the 2019-2020 viral pneumonia COVID-19 outbreak that commenced in Wuhan City, Hubei Province, China (1) (2) (3) (4) . Currently, targeted therapeutics are limited and, at this point in time, effective treatment options remain tentative and highly restricted.

Conventional drug design and drug development strategies require years of investigations and substantial investments. The repurposing of approved pharmaceutical drugs (and drug candidates already in clinical trials) therefore provides an essential, alternative approach to rapidly identify drugs with clinical potential to help manage emerging infectious diseases for which specific treatments and vaccines are not available. A recent high-throughput screening study considered more than 10,000 compounds against the main protease M pro of SARS-CoV-2 (5) . The main protease is one of the key enzymes in the viral life cycle, as it plays a central role in mediating viral replication and transcription. Its structure is shown in Fig. 1 . Along with other, nonstructural proteins (papain-like protease, helicase, and RNA-dependent RNA polymerase) and the spike glycoprotein structural protein, M pro is essential for interactions between the virus and host cell receptors during viral entry (6) , making it an attractive drug target against SARS-CoV-2 (7, 8) .

In a recent study by Jin et al. (5) , a wide range of approved drugs, drug candidates, and natural products were screened by an enzymatic inhibition assay to identify a possible M pro inhibitor. Ebselen, an organoselenium compound whose structure is shown in Fig. 1 , emerged from that study as a promising drug lead to target this crucial enzyme. It was found to exhibit strong antiviral activity in cellbased assays (5) . In previous, pre-COVID-19 work, ebselen was found to have anti-inflammatory, antioxidative, and cytoprotective properties. It has been investigated in the context of multiple diseases, such as bipolar disorders (9) and hearing loss (10, 11) . In addition, a recent report indicated that ebselen also showed potent bactericidal activity against multidrug-resistant (MDR) clinical isolates of Staphylococcus aureus (12) . Likewise, the use of silver and ebselen in synergistic formulations has been shown to be effective against five, clinically difficult-to-treat MDR Gram-negative bacteria (13) . Ebselen has low cytotoxicity (median lethal dose in rats >4600 mg/kg, per os) (14) , and its safety in humans has been evaluated in several clinical trials (10, 11, 15) . Together, these reports-and the underlying data-underscore the clinical potential of ebselen for CoV-2 treatment.

In an effort to understand the molecular mechanisms through which ebselen interacts with M pro , in this work, we present results of atomistic molecular simulations that provide useful, previously unknown insights into the M pro -ebselen complex, which might provide novel avenues to rationally enhance ebselen's activity. First, we present an analysis of the most probable interaction sites between M pro -ebselen, as well as absolute binding free energy calculation for a complementary quantitative evaluation. Second, we evaluate the role of different binding sites on molecular stiffness and molecular strain. Last, we examine how ebselen binding modifies the structure and transport of water in M pro 's catalytic site. We find that ebselen binds to two sites, as opposed to only one at the catalytic site, and we also identify an allosteric mechanism that influences the catalytic site when ebselen is bound at the distant site. Our findings are summarized in a concluding section, where we also provide several suggestions for future experimental work.

Molecular dynamics (MD) simulations of M pro and ebselen were carried out using the AMBER18 simulation package (see Materials and Methods for details). A total of 3 s of atomistic MD trajectories of M pro using ebselen as a molecular probe were organized into 15 replicas of 200 ns each. These trajectories were analyzed to construct probability density maps for the preferred locations of ebselen around the protein. The results are shown in Fig. 1 . Two distinct, highly probable binding sites emerged from this analysis. The first is located within the catalytic site, and the second is in a region that is essential for M pro 's dimerization (16) , between the II and III domains. 1 Similar observations were reported in a previous study (17) , where simple organic solvents or compounds such as acetonitrile, benzene, dimethyl sulfoxide, methanol, phenol, and urea were used as molecular probes. In addition, the authors considered the potential mutability of residues belonging to the M pro catalytic site and explained that the development of drug resistance associated with the natural evolution of M pro could wipe out efforts that target this protein for COVID-19 treatment. Instead, they emphasized the promise of alternative strategies aimed at targeting the region between the II and III domains, which is implicated in dimer formation. The results shown in Fig. 1 are therefore encouraging in that ebselen appears to target both the catalytic site and the dimerization domain. Figure 2 provides a closer look at ebselen's binding motifs in both regions (the catalytic site and the intersection between domains II and III). In the catalytic site ( Fig. 2A) , hydrogen bonding interactions are formed between the carbonyl oxygen of ebselen and the Asn 142 and Gln 189 side chains. One can also appreciate the hydrophobic contacts between ebselen and Met 165 , Pro 168 , Met 49 , His 164 , and both residues belonging to the catalytic dyad, His 41 and Cys 145 . Figure 2B shows a representative configuration of ebselen at the intersection between domains II and III. There are distinct hydrophobic contacts with Phe 294 , Pro 108 , Ile 200 , Val 202 , His 246 , Thr 292 , Ile 249 , Pro 132 , and Ile 249 (some residues are not shown for clarity). Highly dynamic hydrogen bonds were observed between ebselen and the Gln 107 , Gln 110 , and His 246 side chains.

The results shown in Fig. 1 were generated on the basis of direct MD simulations. To arrive at a quantitative estimate of the binding affinity of each site, we used thermodynamic integration (TI) to determine the absolute binding free energy for the catalytic site and the domains II and III site. More specifically, we used particle mesh Ewald molecular dynamics (PMEMD) as implemented in AMBER18 (18) with 11 windows per integration and 10 ns per window. In addition, multiple runs starting from the most probable binding cluster identified in the 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 Fig. 1 , the absolute binding free energies corresponding to both sites are negative, serving to underscore the thermodynamic stability of the ebselen-M pro 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 than that at the catalytic site. A previous report suggested that this small molecule could also inhibit M pro through noncovalent binding, particularly because it has been found to exhibit a stronger inhibition effect than other compounds that were also able to quantitatively modify Cys 145 in the catalytic dyad (5). The results presented here could explain ebselen's high enzymatic inhibition, even when this compound could only partially modify Cys 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 nontrivial, functional conformational changes from nonfunctional 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. It is also important to highlight that different measures of strain have recently been used in studies of allostery in proteins (19) (20) (21) . Here, we used a method originally introduced by our group in the context of local strain in polymeric glasses (22) . Figure 3 shows the shear strain throughout the protein upon binding of ebselen at different sites, alongside the  factor (estimated from the root mean squared fluctuation), which measures thermal fluctuations for each residue. The shear strains are measured relative to the average conformation from all frames of the apo protein trajectory; these strains, therefore, correspond to the deformation upon binding relative to the apo protein. The shear strain for the apo structure is then showing an average internal strain in the protein.

From the  factor analysis, it is evident that when ebselen is bound to the domains II and III interface (Fig. 3A , blue lines), the dynamics of the 44-52 loop, which flanks the catalytic site, is notably 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 to 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 M pro -apo protein (Fig. 3A , green line). The strain analysis helps disentangle the effects of functional and nonfunctional fluctuations and provides a more detailed view of the effects of ebselen binding. As seen in Fig. 3A for the  factor analysis, when ebselen is bound to the domain interface, it produces a large strain at the 44-52 loop (Fig. 3B , blue line). Likewise, the 185-201 loop that also flanks the catalytic site, as well as the region comprising residues 137 to 140, exhibits a high strain, which is not immediately apparent in the  factor. When both molecules are bound simultaneously to M pro (Fig. 3B , black line), a high strain signal is also exhibited around residues 137 to 140, as expected on the basis of the aforementioned  factor results. Last, when ebselen is located in the catalytic site, a much lower strain is showed around all three regions (Fig. 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 to 25). A highly strained region also appears around residues 137 to 140. To understand why the 137 to 140 residues show such a large strain when ebselen is bound to both sites, we turn our attention to the molecular images shown in Fig. 4 .

In the close-up of the residues shown in Fig. 4A , a specific backbone hydrogen bond forms between Lys 137 and Phe 140 because of the binding of ebselen between domains II and III; this conformational change induced by the presence of ebselen causes the high strain shown in Fig. 3B . This conformational change takes place in the middle region between the catalytic site (Cys 145 ) and the binding cleft between domains II and III (Pro 132 ), which points to the relevance of these residues.

From these results, it is evident that when ebselen binds between the II and 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 to 140, where a specific backbone hydrogen bond forms between Lys 137 and Phe 140 . The exact role of this conformational change will be the subject of future studies, but the results presented here show that it serves as a relay between domain III and the catalytic site.

Given that conformational changes in the catalytic site are observed when ebselen binds to M pro far away from this specific region, we turn our attention to the hydration characteristics of the catalytic site in the M pro -ebselen complex (when bound to the domains II and III interface) and in the M pro -apo structure [Protein Data Bank (PDB) code: 6M03] (23). AQUA-DUCT 1.0.5 (24) was used to analyze the water structure and water flux in the M pro protein, with a time window of 50 ns and sampling every 1 ps. The results are shown in Fig. 5 . Figure 5 (top) shows that ebselen binding between domains II and III leads to fewer water inlets compared with the M pro -apo state; this is indicative of a water flux reduction in the catalytic site upon ebselen binding, although it binds far from the catalytic site. From the lower panel in Fig. 5 , it is evident that a volume reduction of the catalytic site occurs when ebselen is located between domains II and 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 molecule that forms a catalytic triad together with Cys 145 and His 41 (17, 25) , and this catalytic water (red in Fig. 5A, bottom) is preserved and remains close to His 41 in the simulations. In contrast, it is clear from Fig. 5B (bottom) that the presence of ebselen induces a displacement and reconfiguration of water hotspots, including the catalytic water (red). These effects could prevent the normal enzymatic function of M pro , as this catalytic water displacement might damage the catalytic triad required for protein activity (26) . Similarly, the observed pocket volume reduction might affect the accessibility of the polyprotein that M pro cleaves, thereby reducing enzymatic function.

As mentioned before, M pro is an attractive drug target against the COVID-19 virus due to its central role in the viral life cycle. A previous structural and evolutionary investigation suggested that the SARS-CoV-2 M pro is not a suitable target for de novo development of inhibitors or the repurposing of drugs against the previous SARS coronavirus (17) . However, only the flexibility and plasticity of the active sites in M pro for COVID-19 and the highly similar previous SARS-CoV M pro were compared. Nevertheless, major differences in both shape and size of the catalytic site were observed, indicating that repurposing SARS drugs for COVID-19 may not be effective.

On the basis of their evolutionary analysis, these authors also pointed out that the virus's mutability will pose further challenges to treatments against the COVID-19 M pro 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 M pro and ebselen. One is located within the catalytic cavity and, more importantly, in the region between the II and III domains, which is essential for M pro dimerization (16) . Detailed calculations of the free energy reveal a higher binding affinity of ebselen to domains II and III than to the catalytic site.

The strain analysis reveals that ebselen bound between the II and III domains exert a pronounced allosteric effect that affects the loops regulating access to the catalytic site. In addition, it also affects residues 137 to 140, where a specific backbone hydrogen bond forms between Lys 137 and Phe 140 .

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 M pro protein.

Together, the results from our comprehensive study of ebselen activity and the underlying data serve to underscore the potential of ebselen for COVID-19 treatment. The discovery of a distant binding site is also encouraging in that it offers potential for development of novel M pro inhibitors based on the ebselen scaffold; indeed, distant sites could be effective targets for other drugs that might have been deemed futile for catalytic sites. Note that current, massive virtual screening campaigns are focused primarily on targeting the SARS-CoV-2 M pro catalytic site; our study, however, highlights the potential importance of a different, distant site (27) . Further experimental structural characterization will be necessary to validate the predictions presented in this work, particularly those pertaining to the distant binding site. In the next stages of this project, we will seek to investigate the experimental crystal structure of the SARS-CoV-2 M pro -ebselen complex, to validate our theoretical predictions and help elucidate the inhibitory mechanism. These efforts could also be combined with rational mutagenesis of interacting residues to alanine and binding assays such as surface plasmon resonance, with which to quantify binding and provide a basis for comparisons to theoretical predictions.

A previous evolutionary study (17) showed that mutating a few residues belonging to the catalytic site is energetically unfavorable. Residues such as P39, R40, P52, G143, G146, or L167 could therefore be considered as key anchoring sites for future M pro inhibitor design. A recent study (5) identified ebselen as a promising repurposed drug against SARS-CoV-2. Our results build on that work by providing a detailed molecular picture of the interactions between ebselen and M pro . The insights provided here could help facilitate the rational molecular design of new analogs, based on the ebselen scaffold and its binding to a distant site, for future coronaviruses. This could be substantial, particularly given the fact that the SARS-CoV-2 M pro protein shares a high sequence similarity to the SARS-CoV M pro (28) , and is therefore a likely target for future coronaviruses. 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 molecules. We conclude by pointing out that our focus in this work has been on noncovalent complexes between ebselen-M pro . In a subsequent stage, we also plan to investigate covalent complexes involving Cys 145 .

A total of more than 6 s of classical MD simulations of SARS-CoV-2 M pro -apo state and M pro -ebselen complex were run using the AMBER18 (29) simulation package (3 s, ebselen as a molecular probe; 2.4 s, shear strain analysis; 100 ns, water structure and flux analysis; 990 ns, free energy analysis). The receptor initial configuration for the M pro -ebselen system was taken from a recently reported structure for the M pro -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 (30, 31) . The partial atomic charges were determined by the restrained electrostatic potential fitting technique. Those electrostatic potential calculations were performed at the HF/6-31G level with Gaussian 09. For M pro -apo simulations, we used the recently reported structure (23) (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 (32) . 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/mol·Å 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 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 nonbonded cutoff of 10 Å, and the SHAKE algorithm was used to implement rigid constraints. 

To determine the most probable interaction sites between M pro -ebselen (global hotspot), 15 replicas of 200 ns each were run (a total of 3 s), where configurations were saved every 20 ps. The CPPTRAJ (33) software was used to process the trajectories, where at the first stage all trajectories were aligned by means of minimizing the distance among protein backbone atoms (C, N, and CA). Grid command was used to track the ebselen molecule and to produce a number density map, where the grid resolution was selected to be 0.5 Å.

Initial M pro -ebselen complex configurations were selected from previous global binding affinity study, with ebselen being bound to the catalytic site, to the intersection between domains II and 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 M pro -apo structure (PDB code: 6M03), and three 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 (19-21, 34, 35) . A radius R around each central atom i containing n other atoms j defines the local neighborhood around the central atom. The instantaneous position of atom i at any timestep in the MD is x i , and the position of the same atom at any timestep of the reference simulation is x 0, i . To first order, the distances between atom i and its neighbors j are related through the deformation matrix F by x j − x i = F(x 0, j − x 0, i ), which forms an overdetermined system of linear equations, and an optimized F* is sought by minimizing the difference between the actual distances and the projected distances to an affine deformation: min  j=1 n [ x j − x i − F( x 0,j − x 0,i ) ] 2 . The atomic strain tensor is then found by  = 1 _ 2 ( F T F − I) , and the magnitude is defined as the L2 norm of the shear part of the strain tensor: Tr(  2 ) − 1 _ 3 (Tr( )) 2 , since proteins are generally incompressible (19, 21) . For this analysis, a radius of 10 Å around each atom is considered (21) , 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 M pro protein of SARS-CoV-2, we used AQUA-DUCT 1.0.5 (24) . We obtained inlets, MAV, and hotspots for water molecules in the M pro 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, and D187) (17) . The M pro protein was studied in two different scenarios, with ebselen in the domains II and III site and the apo protein with no ligand. The time window used in both calculations was 50 ns, sampling every 1 ps. Images were created with open-source PyMOL (36) .

The absolute binding free energy is defined as follows: G binding = G L − G RL , where G RL is the free energy change of ebselen annihilation in the M pro complex and G L is the free energy change of ebselen annihilation in water. To calculate these free energy changes, we use TI implemented in PMEMD for AMBER18. We use the onestep annihilation protocol with soft-core potentials (37) . 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 (λ = 0.1) with 10 ns of simulation time per window. To keep the ligand from wandering in TI calculations, we used a soft restraint of 10 kcal/ mol·Å 2 (37) .

China Novel Coronavirus Investigating and Research Team, A novel coronavirus from patients with pneumonia in China

Early transmission dynamics in wuhan, china

A new coronavirus associated with human respiratory disease in China

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

Coronaviruses -drug discovery and therapeutic options

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 noise-induced hearing loss

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

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

Structural and evolutionary analysis indicate that the SARS-CoV-2 Mpro is an inconvenient target for small-molecule inhibitors design

Improving the efficiency of free energy calculations in the amber molecular dynamics package

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

A deformation gradient tensor and strain tensors for atomistic simulations

Proteins: The physics of amorphous evolving matter

Calculation of local mechanical properties of filled polymers

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

Coronavirus main proteinase (3CL pro ) 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

Research and development on therapeutic agents and vaccines for COVID-19 and related human coronavirus diseases

Automatic atom type and bond type perception in molecular mechanical calculations

Development and testing of a general amber force field

Simmerling, ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB

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

Pymol: An open-source molecular graphics tool. CCP4 Newsletter on protein crystallography

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

We are grateful to the Research Computing Center of the University of