key: cord-0974358-z09dqal2
authors: Santra, Santanu; Giri, Santanab; Jana, Madhurima
title: Unraveling the Origin of Interactions of Hydroxychloroquine with the Receptor-Binding Domain of SARS-CoV-2 in Aqueous Medium
date: 2020-12-19
journal: Chem Phys Lett
DOI: 10.1016/j.cplett.2020.138280
sha: 6c993b64ab15e03d91984c6cc84f7af78a60178d
doc_id: 974358
cord_uid: z09dqal2

Interactions of hydroxychloroquin (HCQ) with the receptor binding domain (RBD) of SARS-CoV-2 were studied from atomistic simulation and ONIOM techniques. The key-residues of RBD responsible for the human transmission are recognized to be blocked in a heterogeneous manner with the favorable formation of key-residue:HCQ (1:1) complex. Such heterogeneity in binding was identified to be governed by the differential life-time of the hydrogen bonded water network anchoring HCQ and the key-residues. The intermolecular proton transfer facilitates the most favorable Lys417:HCQ complexation. The study demonstrates that off-target bindings of HCQ need to be minimized to efficiently prevent the transmission of SARS-CoV-2.

The devastating nature of novel coronavirus , emerged in December, 2019 causes worldwide millions of rapid deaths and thus on March, 2020 the World Health Organization (WHO) officially declared it as pandemic [1] . The receptor binding domain (RBD) of spike protein of Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pathogen having close genome identity with SARS-CoV was reported to involve in molecular recognition to the host cellular receptor, angiotensin-converting enzyme 2 (ACE2) [2] [3] [4] [5] [6] . Although, a series of rapid diagnostic tools are available for the fast detection of the disease, to the date no specific vaccines and therapeutic treatments are available. The development and manufacture of new vaccine against COVID-19 is time-taken whereas searching for therapeutics, like short peptide inhibitors against the virus is relatively faster; however, their applicability is not known [7, 8] . An alternative strategy to find a quick solution to fight against the virus is the clinical trial of novel and repurposed drugs [9, 10] .

Among the antiviral drugs, Chloroquine and Hydroxychloroquine (HCQ) have garnered intense attention and HCQ, a popular drug for the treatment of malaria, rheumatoid arthritis and lupus has been reported to reduce viral load and minimize the SARS-CoV-2 transmission [11, 12] . So far, mixed responses about HCQ against COVID-19 have been noticed and therefore, on 23rd March, 2020 the executive group of WHO implemented temporary pause of HCQ,[13] however, later on 3rd June, 2020, based on the mortality data the committee has recommended its use [14] .

Thus, the situation demands special attention and further investigation on HCQ against the novel coronavirus. In this aspect physical origin of interactions of HCQ with the RBD of SARS-CoV-2 is expected to be crucial. As of now, the mechanism by which HCQ inhibits SARS-CoV-2 is not known. A recent study showed that the structural interactions of RBD of SARS-CoV-2 and ACE2, responsible for human transmission may occur through few specific amino acids (contact residues) present over the extended loop region of RBD [15] . Therefore, it would be of fundamental interest to study the molecular interactions between HCQ and RBD to explore the efficiency of HCQ in the targeted binding of the contact residues of RBD's active-sites for the inhibition. Understanding of such interactions would provide wealth of information on HCQ's effectiveness and limitations to block the contact residues, responsible for the transmission. In this work, the molecular level understanding of interactions and the binding efficiency of the contact residues, hence forth named as key residues; [15] Lys417, Tyr453, (region 1) Gln474, Phe486, (region 2) Gln498, Thr500, and Asn501 (region 3) ( Fig. 1(a) ), of the RBD of SARS-CoV-2 with HCQ in aqueous solution were explored. Efforts have been made to investigate how effectively HCQ binds with the targeted key-residues to inhibit RBD. Additionally, attempts have been made to identify the most efficient key-residue and the physical origin of its binding.

Two different atomistic Molecular Dynamics simulations of the Receptor Binding Domain (RBD) (PDB ID: 6M17) [15] of SARS-CoV-2 were carried out at 300 K temperature in pure water and in presence of HCQ. We have used CHARMM36 all atom force field and potential parameters for the protein [16] and CHARMM general force field for the HCQ [17] . The source of the initial structure and the parameters of HCQ are given in the Supplementary Material. The web based ligand reader and modeler of CHARMM-GUI was used to model HCQ. [18] The CHARMM general force field was found to be used in some of the recent studies for the modeling of HCQ. [19, 20] The modified TIP3P water model compatible with CHARMM36 force field and NAMD code has been used to model water molecules [21] . RBD was immersed in two separate cubic cells of edge length 90 Å containing pure water and aqueous solution of HCQ (~0.01 M). Use of concentrated HCQ solutions around the chosen concentration or higher can be found in several other studies [22, 23, 24] To avoid any unfavorable contacts, the insertion process was carried out by carefully removing all those added solvent molecules that were found within 2 Å from any of the atoms of the RBD. The HCQ molecules were placed randomly in the box. The overall charge of the system was neutralized by adding appropriate number of counter ions. To eliminate initial stress, the systems were first minimized using the conjugate gradient energy minimization method as implemented in the NAMD code [25] . The temperature of the systems was gradually increased to the targeted 300 K within a short MD run under isothermal-isobaric ensemble (NPT) conditions at constant pressure of 1 atm. The systems were equilibrated for 5 ns at the desired temperature under NPT ensemble conditions. Thereafter, long production trajectories at 300 K were generated in NVT ensemble for around 150 ns duration. All the simulations were carried out using NAMD code with a time step of 1 fs. Temperature of the systems was controlled by using Langevin dynamics method with a friction constant 1 ps −1 . The pressure of the systems was controlled by Nosé-Hoover Langevin piston method [26] . The minimum image convention [27] was employed to calculate the short range Lennard Jones interactions using a spherical cut-off distance of 12 Å with a switch distance of 10 Å. The long range electrostatic interactions were calculated using the Particle-Mesh-Ewald (PME) method [28] .

The electronic structure of binding interaction of HCQ with the key-residue of RBD was quantified by computing the binding energy of key-residue:HCQ (1:1) complex by performing a two-layer ONIOM (QM:MM) [29] optimization in Gaussian09w [30] . Initially, HCQ was placed within the distance of 3.5 Å from each of the key-amino acid residue of RBD. This is essentially within the hydrogen bond distance as formation of hydrogen bond between drug-receptor plays crucial role in the inhibition process. In this framework universal force field [31] was used for optimization of MM level (lower layer) that consists of entire RBD except the particular key residue that forms complex with HCQ, whereas for the optimization of QM part (higher layer), consists of key amino acid residue and HCQ, wB97XD [32] hybrid functional with 6-31G(d,p) basis set was used. The binding energy was calculated by using the following equation:

Where is the energy of total complex, is the energy of the RBD computed from a two-layer ONIOM calculation where key residue was treated in high layer (QM), and E HCQ is the energy of the HCQ molecule calculated at wB97XD/6-31G(d,p) level of theory. A post Hartree-Fock method, MP2 [33] was additionally used to calculate binding energy on wB97XD/6-31G(d,p):UFF optimized geometry. As water plays important role in biological processes, binding energy was also calculated by employing Polarizable continuum model (PCM) [34] . To reduce computational cost arising from the large size of RBD, the solvent phase binding energy calculations were performed by considering only key-residue-HCQ (1:1) complex without changing its geometry in the protein environment.

Considering the crucial role of protein's conformational flexibility in modulating drug binding kinetics, the flexibility of SARS-CoV-2's RBD conformation was quantified by computing the root-mean-square-deviations (RMSD), and the residue wise root-meansquare-fluctuations (RMSF) in the experimental HCQ solution, and the results were compared with that obtained from the RBD in pure water simulation. Apparently, from Fig.   1 (b) the flexibility of the RBD is observed to be nearly similar in pure water as well as in HCQ solution; however quantification of residue wise RMSF (Fig. 1(c) ) shows that apart from few residues of binding region 2, most of the residues including the key-residues fluctuate less in presence of HCQ. Such relative restriction of the RBD's residue false-positive results that may influence off-target inhibition and raises toxicity [35] . To avoid such common problem in drug discovery, it is necessary to diminish the off-target interactions of HCQ (logP ~ 4 and topological polar surface area, TPSA ~ 48.38 [36] ) by controlling their selfaggregation behaviour. Development of simple derivative of HCQ with lesser logP and higher TPSA (<80) might be helpful to prohibit the self-association of the drug like HCQ [37, 38] . As the principal interactions generally appear from the sites of the key-residues, we have quantified the probability of at least only one HCQ molecule to be present within the solvation shell of the individual key-residue of the RBD over the equilibrated trajectory. To identify 1:1 keyresidue:HCQ complex, we imposed a distance cut-off criteria. According to the criteria, when a HCQ molecule was found to be present within the distance of 5 Å, from each of the key residue then that key-residue and the tagged HCQ was considered to form 1:1 complex. In Figure S -III of the Supplementary Material we have shown the representative 1:1 complexes extracted from the equilibrated trajectory. As depicted in Table I , among all, Lys417, and Tyr453, i.e the key-residues from region 1 have high probability to form 1:1, (key-residue:HCQ) complex over the simulation period. The lesser probability of forming such complex by the other key-residues with HCQ infers that although HCQ remains within the solvation shell of the binding region (discussed above) of RBD, they do not form complexes with the targeted sites. In this context, it is noteworthy that there exists probability of formation of 1:2 complexes, however as the probability was observed to be very less (~ 0.06% -4%),

we have not considered 1:2 complexes for the further investigation.

Next, to explore the change in RBD-water hydrogen bonding network due to the presence of HCQ in the RBD's hydration shell we have measured the life-time of the hydrogen bonds formed between the key-residues and water and compared their life-time when no HCQ is present in the solution. A pure geometric criterion was used to define hydrogen bonds [39, 40] . According to these criteria the first condition for an atom of the protein to form a hydrogen bond with the solvent is that the distance between the tagged atom and the oxygen atom of the water molecule with which it is hydrogen bonded be of such key-residue-water hydrogen bonds are also shown in the Table I . Presence of longlived hydrogen bonded water of this kind clearly indicates that the key-residues and HCQ are anchored by this type of water molecules (Fig. 2 (d) ), and such trapped water generally survives for a long time due to the confinement. This makes the drug binding process feasible at the RBD surface. At this point, it is important to quantify the binding interactions between the key-residues and HCQ to explore whether the identified 1:1 interaction is energetically favourable or not and among all, which keyresidue could play significant role. Thus, we studied the electronic structure of the binding site HCQ:key-residue complex by adopting ONIOM (QM:MM) methodology to compute binding energy of key-residue:HCQ (1:1) complex in the protein environment. Fig. 3 shows the optimized geometries of the pairs. The calculated binding energy as obtained from the two different levels of theories (Table II) in gas phase reveals that except Phe486, all the key-residues bind quite efficiently and favourably with HCQ. 

In summary, HCQ was noticed to bind favourably but heterogeneously with all the targeted key-residues of RBD of SARS-CoV-2 by forming 1:1 complex. The most favoured binding between HCQ and Lys417 was found to be governed by the intermolecular proton transfer phenomenon. The intercalated long-lived hydrogen bonded water molecules present in the shared solvation regions of HCQ and key-residues were identified to play important role in key-residue:HCQ complexation process. It could be relevant to study the local properties by using Density Functional Theory treatment that may provide more signals in the modifications of the binding energies of key-residue-HCQ. Further, the HCQ solution at various concentrations is likely to provide wealth of information and physical insights about the drug's potency to bind SARS-CoV2. Work in both the directions is ongoing.

However, we believe that the mechanistic knowledge and the information obtained from the study in general, would facilitate the development of effective potential therapeutics for the prevention of activation of SARS-CoV-2 virus.

A Pneumonia Outbreak Associated with a New Coronavirus of Probable Bat Origin

Structure of SARS Coronavirus Spike Receptor-Binding Domain Complexed with Receptor

Cryo-EM Structure of The SARS Coronavirus Spike Glycoprotein in Complex with Its Host Cell Receptor ACE2

Cryo-EM Structure of the 2019-nCoV Spike in the Prefusion Conformation

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

Computational Design of ACE2-Based Peptide Inhibitors of SARS-CoV-2

Molecular Modeling and Chemical Modification for Finding Peptide Inhibitor Against Severe Acute Respiratory Syndrome Coronavirus Main Proteinase

Living Systematic Review Protocol for Covid-19 Clinical Trial Registrations

World Health Organisation: Coronavirus disease (COVID-2019) situation reports

In Vitro Antiviral Activity and Projection of Optimized Dosing Design of Hydroxychloroquine for the Treatment of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2)

Hydroxychloroquine, A Less Toxic Derivative of

Is Effective in Inhibiting SARS-Cov-2 Infection in Vitro

Solidarity" clinical trial for COVID-19 treatments

Structural Basis for the Recognition of the SARS-CoV-2 by Full-Length Human ACE2

Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles

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

CHARMM-GUI Ligand Reader and Modeler for CHARMM Force Field Generation of Small Molecules

Inhibitory activity hydroxychloroquine on COVID-19 main protease: An insight from MD-simulation studies

Synergistic antiviral effect of hydroxychloroquine and azithromycin in combination against SARS-CoV-2: What molecular dynamics studies of virus-host interactions reveal

Solvent-Induced Forces between Two Hydrophilic Groups

A Novel RP-HPLC-DAD Method Development for Anti-Malarial and COVID-19 Hydroxy Chloroquine Sulfate Tablets and Profiling of In-Vitro Dissolution in Multimedia

Inhaled hydroxychloroquine to improve efficacy and reduce harm in the treatment of COVID-19

Scalable molecular dynamics with NAMD

Constant pressure molecular dynamics simulation: The Langevin piston method

Computer Simulations of Liquids

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

A New ONIOM Implementation in Gaussian98. Part I. The Calculation of Energies, Gradients, Vibrational Frequencies and Electric Field Derivatives

Gaussian 09, Revision A.02

UFF, A Full Periodic-Table Force-Field for Molecular Mechanics and Molecular-Dynamics simulations

Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections

Direct MP2 Gradient Method

Quantum Mechanical Continuum Solvation Models

Compound Aggregation in Drug Discovery: Implementing a Practical NMR Assay for Medicinal Chemists

Fast Calculation of Molecular Polar Surface Area as A Sum of Fragment Based Contributions and Its Application to The Prediction Of Drug Transport Properties

Comparative Molecular Dynamics Simulation of Aggregating and Non-Aggregating Inhibitor Solutions: Understanding the Molecular Basis of Promiscuity

Computational Analysis of Calculated Physicochemical and ADMET Properties of Protein-Protein Interaction Inhibitors

Molecular Dynamics of Water at the Protein Solvent Interface

Residue Specific Interaction of an Unfolded Protein with Solvents in Mixed Water−Ethanol Solutions: A Combined Molecular Dynamics and ONIOM Study

Ligand-induced Proton Transfer and Low-Barrier Hydrogen Bond Revealed by X-ray Crystallography

Spontaneous Proton Transfer to a Conserved Intein Residue Determines on-Pathway Protein Splicing

The work carried out by using the computational facility created under the grant no.EMR/2017/001325, Science and Engineering Research Board (SERB) Government of India.