Unraveling the origin of interactions of hydroxychloroquine with the receptor-binding domain of SARS-CoV-2 in aqueous medium.

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.

Introduction

The devastating nature of novel coronavirus (COVID-19), 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.

Fig. 1

(a) RBD (yellow) of SARS-CoV-2 along with all the binding site residues (red). The key-residues present in three different regions are shown separately. (b) Time evolution of RMSDs of the whole RBD (heavy atoms) in pure water and in HCQ solution. (c) Residue-wise RMSF values of the RBD. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Computational details

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 EHCQ 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.

Results and discussion

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-mean-square-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 wise fluctuation infers that in presence of HCQ, reorganization of the solvation shell of the RBD occurs that hinders the fluctuation of the residues of the RBD.

Solvation shell reorganization is extremely important for biomolecular recognition and protein-drug binding. As the key-residues of extended loop region of RBD were reported to directly participate in H-bonding interactions with the binding domain of ACE2 [7], [15], reorganization of the solvation layer in presence of HCQ is essential for the establishment of RBD-HCQ direct hydrogen bonding network or for the feasible formation of water-mediated hydrogen bonding network to inhibit the anchoring sites of the RBD, through which it generally interacts with the host protein. We observed that the average number of HCQ, (quantified from the equilibrated trajectory) present in the solvation shell of the binding region residues of the RBD (~8), as compared to that around whole RBD (~14) is more than 50%. This primarily infers that the HCQ molecules remain close to the residues of binding regions of the RBD and hence can take part in the complex formation process. The snapshot of RBD and its surrounded HCQ (Fig. 2 (a) and (b)) as obtained from the simulated configuration along with its starting configuration clearly depicts this fact. However, it can be importantly noticed from the figure that the HCQ molecules self-aggregate on the RBD surface. It may be further noted that the self-aggregation phenomenon of HCQ in the solution at the studied concentration was additionally observed from the separate simulation (~100 ns) performed for HCQ in aqueous solution in absence of the RBD. The snapshots of the simulated configurations as shown in Figure S-I of the Supplementary Material indicate HCQ’s aggregation phenomenon in solution. Such aggregation phenomenon of HCQ could be a concentration dependent phenomenon; in-depth studies need to be performed to address this issue. However, our study reveals that HCQ can self-aggregate to form clusters of different sizes involving 2–5 members. A representative snapshot of the HCQ cluster near the RBD surface is also shown in Fig. 2 (c). Additionally, the arrangements of HCQs around the RBD at the regular time interval as obtained from the equilibrated trajectory are shown in the Figure S-II of the Supplementary Material. Such self-aggregation of drugs might provide 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 self-aggregation 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 key-residue: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 1 , 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.

Fig. 2

(a) The initial position of the HCQ molecules in solution and (b) the position of HCQ at the end of the simulation. Water molecules are not shown for clarity. (c) A representative snapshot of self-aggregated HCQ linked through hydrogen bonds. (d) A representative water bridged hydrogen bonded network formed between Lys417 and one of the HCQ present in the solution.

Table 1

Probability of forming key-residue:HCQ, (1:1) complex over the simulation period and the key-residue-water hydrogen bond life times as obtained from the equilibrated trajectory.

Residue	% Probability of forming key-residue:HCQ (1:1) complex	τavg (ps)		τavgshared (ps)
		Water	HCQ solution
Lys417	91.4	3.98	6.08	49.10
Tyr453	87.8	3.09	5.12	20.33
Gln474	6.13	2.56	2.81	17.47
Phe486	27.4	2.85	4.0	30.13
Gln498	10.8	1.95	2.33	19.76
Thr500	8.6	2.40	3.08	22.03
Asn501	6.4	2.45	2.83	20.78

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 within 3.5 Å. The second condition is the donor-hydrogen-acceptor angle should be within 140°. The average life-time (τavg) values as shown in Table 1, clearly infer that the presence of HCQ increases the life-time of the key-residue-water hydrogen bonds. Among all, longer life-time of Lys417-water hydrogen bond is prominent from the table. The intercalated water molecules present within the shared solvation region of a protein and ligand are specifically crucial in biomolecular recognition since they can form bridged hydrogen bonds between them. Therefore, the life-time of the key-residue-water hydrogen bonds were additionally computed by considering the water molecules that were present only at the shared solvation region of the key-residue and HCQ. The average life times (τavg shared) of such key-residue-water hydrogen bonds are also shown in the Table 1. Presence of long-lived 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 key-residue 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 2 ) in gas phase reveals that except Phe486, all the key-residues bind quite efficiently and favourably with HCQ.

Fig. 3

Ground state optimized geometries of HCQ bound with the key-amino acid residues (1:1 complex) of RBD of SARS-CoV-2.

Table 2

Binding energies (kcal/mol) of key-residue-HCQ, 1:1 complex obtained by using ONIOM(QM:MM) two-layer technique.

Complex	Eb (Gas)	Eb (Gas)	Eb (H2O)	Eb (H2O)
HCQ@	wB97XD/6-31G(d,p):UFF	MP2/6-31G(d,p):UFF	wB97XD/6-31G(d,p)	MP2/6-31G(d,p)
Lys417	−113.17	−80.86	−93.14	−89.59
Tyr453	−124.81	−67.11	−6.53	−6.48
Gln474	−45.53	−82.21	−5.11	−4.62
Phe486	7.69	18.47	−11.15	−9.45
Gln498	−44.51	−69.86	−7.05	−6.77
Thr500	−50.78	−74.93	−8.84	−7.30
Asn501	−45.05	−63.66	−7.94	−7.00

It may be noted that the reported binding energies are the most favourable one, relative to the binding energies obtained from the several trials of other optimized geometries of the complexes constructed individually. Additionally, the binding energies were computed for several identified 1:1 complexes those were extracted from the equilibrated trajectory of MD simulation. The average binding energies obtained for each of such complexes are given in the Table S-IV of the Supplementary Material. Importantly, we observed that, the 1:1 complexes identified from either of the techniques possess favourable binding interactions. Since water plays an important role in protein-drug interactions, the binding energies were additionally calculated in water medium using PCM model [34]. However, to reduce computational cost, the interactions were calculated considering only the key-residue and HCQ, without altering their optimized structure in the protein environment. It is interesting to observe from Table 2 that the unfavourable binding of Phe486 with HCQ in gas phase now becomes favourable in presence of water medium, implying crucial role of water in protein-drug binding in general. The binding in aqueous environment was also found to be favourable for all the other pairs. Remarkably, among all, Lys417 was observed to possess higher binding affinity towards HCQ in aqueous phase, too. Consistency of Lys417′s higher binding affinity towards HCQ, further inspired us to look for the possible reason of such affinity. The closer view of the optimized geometry of Lys417-HCQ reveals that there is a possibility of transfer of H between the two N centers of Lys417 and HCQ, as one of the H of Lys417-NH3 is observed to be detached from its parent atom and initiates to form a bond (1.06 Å) with the pyridinium (Py) N of HCQ. This is an important observation, and could be the origin of the high binding affinity of Lys417 towards HCQ.

To observe the feasibility of the Py@N-H bond formation, we have performed an energy scan of N-H bond of lysine keeping the rest of the geometry intact. Fig. 4 depicts the change of energy along the N--H distance in gas as well as in aqueous phase. It clearly demonstrates that Py@N-H bond formation is a favourable process (reaction energy, ΔER (Gas) = −5.52 kcal/mol and ΔER (H2O) = −3.68 kcal/mol) with a very narrow activation energy barrier (ΔE#) of 0.7 kcal/mol in gas phase and 2.25 kcal/mol in aqueous phase. Such low energy barrier or barrier less proton transfer phenomenon in protein environment has been observed in several other studies [41], [42]. A movie that illustrates proton transfer has been shown in the Supplementary Material. To ensure the possible transition state (TS) structure, the conformation having maximum energy in the energy scan diagram was extracted for consequent vibrational frequency analysis. The existence of one imaginary frequency associated with the H translation between two N atoms of Lys-NH3 and Py-N of HCQ confirms the fact further.

Fig. 4

A possible H-transfer reaction between Lys417 and HCQ.

Conclusion

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.

CRediT authorship contribution statement

Santanu Santra: Formal analysis. Santanab Giri: Methodology, Investigation, Formal analysis, Writing - original draft. Madhurima Jana: Conceptualization, Methodology, Formal analysis, Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.