Encapsulated hydroxychloroquine and chloroquine into cyclic oligosaccharides are the potential therapeutics for COVID-19: insights from first-principles calculations.

Novel-Coronavirus (COVID-19) outburst has become a worldwide pandemic which threaten the scientific community to design and discover efficient and effective treatment strategies against this deadly virus (SARS-CoV-2). Still now, there is no antiviral therapy or drug available in the market which can efficiently combat the infection caused by this virus. In this respect, using available drugs by screening with molecular docking and molecular dynamics studies not only minimizes lengthy chemical trials but also reduces discovery cost for the pharmaceutical industry. During the COVID-19 pandemic situations hydroxychloroquine, chloroquine known as HCQ and CQ tablets have gained popularity as for the treatment coronavirus (COVID-19) but the main threatening effect of HCQ, CQ use lies on their side effects like blistering, peeling, loosening of the skin, blurred vision stomach pain, diarrhea, chest discomfort, pain, or tightness, cough or hoarseness which require immediate medical attention. Encapsulation of HCQ and CQ drugs by the cyclic macromolecules such as α and β-Cyclodextrin, to form host-guest complexes is very effective strategy to mask the cytotoxicity of certain drugs and alleviating and modulating side effects of drug applications. In the present work, we have encapsulated the HCQ and CQ drugs α and β-Cyclodextrin and made a comprehensive analysis of stability, optical properties. Details analysis verified that between QC and HCQ, HQC showed stronger affinity towards β-Cyclodextrin. This strategy can reduce the side effect of HCQ and CQ thereby offers a new way to use these drugs. We hope the present study should help the researchers to develop potential therapeutics against the novel coronavirus.

Introduction

In late December 2019, a lethal virus viz. SARS-CoV-2 has unfold quickly worldwide making a demoralizing impact on hundreds of millions of lives of approximately 210 countries [1], [2], [3]. Due to its high mortality rate WHO confirmed a pandemic situation after the three months of its preliminary screening [4]. Till now, World Health Organization reported on 13thJuly 2021, the virus had caused 190,923,871 infections and 4027,861deaths all over the world [5]. This pandemic slowdowns the national healthcare systems and global economy. In the lungs, the virus targets cells expressing angiotensin-converting enzyme 2 resulting in the overall decrease in oxygen levels in blood and leaving the patient in fatal condition. COVID-19 shows a spectrum of clinical presentations ranging from asymptomatic to severe respiratory failure with some common symptoms like fever, cough, headache, and diarrhea [6], [7], [8]. Structurally SARS-CoV-2, the single-stranded positive-sense RNA virus is closely connected to the genomic organization of SARS-CoV recognized in 2003 [9]. The virus entry is initiated by the interaction of angiotensin converting enzyme 2 of human cell and the spike protein of SARS-CoV-2 [10]. Immediately after the virus entry “budding”) and viral replication and transcription are started with the functional proteins like main protease (Mpro), RNA-dependent RNA polymerase (RdRp), papain-like protease (PLpro). Though there are four different types of coronaviruses (α-COV, β-COV, γ-COV and δ-COV)are known but the work function of SARS-COV-2 remains unclear [11,12].

Numerous attempts have been made by different groups of scientist either new or repurposed drugs but yet now there is no active drug present in the market which can successfully combat SARS-CoV-2.The leading method for these drug developments include repurposing of approved drugs or the use of different theoretical methodologies like molecular docking, molecular dynamics simulations, virtual screening which could not only save time but cost effective and facilitate the drug discovery of COVID-19.

During the different clinical trials, some well-known drugs exposed intoxicating effect on this disease by curing a few COVID-19 patients. But the actual modes of action of those drugs are still unknown. Recently a number of drugs are reported viz. Remdesivir [7,10], Hydroxychloroquine (HCQ) [11], Chloroquine (CQ) [13], and natural products of different pant resources [14], [15], [16], [17], [18] and other RNA virus drugs [19] but none of these are particularly effective for this virus. Hydroxychloroquine in combination with azithromycine was also applied for the treatment of COVID-19 [20]. HCQ and CQ are lipophilic weak bases which diffuse through cell membranes and into different parts of the cell like endosomes, lysosomes, and golgi vesicles, where they can be easily protonated, and trapped. Currently there are many trials that have been registered globally involving either CQ or HCQ alone or in combination with other drugs such as azithromycin. It was also reported that the addition of azithromycin, an antibiotic, to the hydroxychloroquine arm resulted in complete viral clearance (6/6, 100%) versus hydroxychloroquine alone (8/14, 57%) [21,22]. Chloroquine phosphate has shown positive activity on COVID-19 and also applied to COVID-19 patients [23,24]. Although in COVID-19 pandemic situations hydroxychloroquine, chloroquine known as HCQ and CQ tablets have gained popularity as for the treatment coronavirus (COVID-19), the main threatening side effect blistering, peeling, loosening of the skin, blurred vision stomach pain, diarrhea, chest discomfort, pain, or tightness, cough or hoarseness which require immediate medical attention [25]. For this reason FDA and WHO warns the use of HCQ and CQ and under the present situation there is not enough information to find out other potent drugs with immediate effect. Therefore, a global response with new strategies is urgently needed to effectively use these drugs against this exceptional pandemic disease.

Though there are few reported works, some of which already referred above, most are still in clinical trials and no results from in-vivo experiments have been explored yet. The advancement of encapsulation of different drugs by the cyclic macromolecules such as α and β-Cyclodextrin, to form host-guest complexes is very effective strategy to mask the cytotoxicity of certain drugs and alleviating and modulating side effects of drug applications. In these circumstances the encapsulation strategy can reduce the side effect of HCQ and CQ thereby offers a new way to use these drugs. In the present work, we have encapsulated the HCQ and CQ drugs in α and β-Cyclodextrin and made a comprehensive analysis of stability, optical properties, thermodynamic properties of the host-guest inclusion complexes.

Computational details

All density functional theory (DFT) calculations in the present work were carried out using the Gaussian 16 program [26]. Ground sate geometry optimizations of the of α–CD and β–CD,CQ, HCQ and the inclusion complexes were performed at M06–2X/6–31+G(d)level of theory [27]. Hybrid B3LYP functional reports relatively poor stacking interactions in non-covalently bonded units while the meta-generalized gradient approximation M06–2X, with the 6–31+G(d) basis set are pretty reliable and precisely describes the non-covalently bonded interaction energies (hydrogen bonding, π–π stacking) present in the π-system [28,29]. The optimized geometries correspond to minima on the potential energy surfaces were confirmed by vibration frequency analysis at the same level of theory. Different types of weak interactions like H-bonding, van der Waals, steric interactions were visualized by Non Covalent Interaction (NCI) [30] index plots of the reduced density gradient (RDG or s) vs. molecular density ρ were analyzed using the Multiwfn 2.6 [31] suite at the ground state geometries. Molecular electrostatic potential (MESP) maps were generated at the same level of theory to understand the extent of existing charge transfer interactions in the inclusion complexes. Furthermore, adsorption energies or binding energy (ΔEads) for all inclusion complexes were evaluated by the following formula:

Where EComplex, ECD, EHCQ are the total energy of the geometry optimized HCQ/α-CD or HCQ/β–CD complexes, free α–CD/β–CD and the guest HCQ molecules, respectively. The host guest interaction energies were further rectified using basis set superposition errors (BSSE) employing the counterpoise procedure of Boys and Bernardi [32,33] as implemented in G16 package.

The molecular dynamics of the interaction of CQ/α-CD/β –CD or HCQ/α-CD/β –CD complexes were simulated employing the atom centered density matrix propagation (ADMP) approach at the same level of theory with step-size of 1000 over a timeframe of 300 femtosecond. The geometric parameters and energy fluctuations of these complexes for each 500 fs are analyzed.

Molecular Docking simulations were performed by selecting best binding sites using AutoDockVina package [34]. Auto Dock Tools was used to construct the starting structure.

Results and discussion

Conformational analysis of the host-guest complexes

Fig. 1 illustrates the graphical demonstration of the energy changes associated throughout the inclusion process of HCQ and CQ drugs in α and β-Cyclodextrin cavities at different Z positions for two different orientations.

Fig. 1

Binding energies profile of (a) CQ-α–CD (b) CQ-β-CD (c) HCQ-α–CD (d) HCQ-β-CDcomplex at different positions Z accounting basis set superposition error (BSSE).

The guest is moved through the host cavity along the Z-axis from −4.0 to +10.0 Å with a stepwise of 1.0 Å from right side as well as left side as indicated in Scheme 1 . At each step the binding energies of the host-guest complexes are computed and plotted in the Fig. 1 by orientation A and B with black and red colors respectively.

Scheme 1

Orientation A and Orientation B for CQ-α-CD, CQ-β–CD, HCQ-α–CD and HCQ-β-CD inclusion complexes.

A detailed analysis of the Fig. 2 revealed that the orientation A always yields greater binding energies than the orientation B indicating that the host-guest inclusion complexes are preferentially formed by the orientation A fashion. Our simulated results agreed well with the experimental results obtained from Roy et al. [35].

Fig. 2

Geometry optimization of the molecules at M06–2X-/6–31+G(d) level of theory. Top and side views of the (a) CQ-α–CD (b) CQ-β-CD (c) HCQ-α–CD (d) HCQ-β-CD inclusion complexes respectively.

Stable geometries and adsorption energy analysis

Ground state optimized geometries of the studied host-guest inclusion complexes are illustrated in Fig. 2. The mode of binding of CQ and HCQ with α–CD/β–CD is distinctly different. Due to the presence of additional stronger H-bonding of HCQ with CD ring hydroxyl group, the HCQ moiety is not deeply penetrated within the cavity of α–CD/β –CD hosts (Fig. 2) yielding greater binding energies with respect to corresponding CQ complexes. Due to the bigger cavity size of β-CD (7.80 Å) over α-CD (5.70 Å) [36] CQ and HCQ moieties showed deeper penetration into the β-CD cavity.The calculated binding energies along with global and geometric parameters like HOMO, LUMO, bang gap, global hardness, global softness, electro-negativities are illustrated in Table 1 in vacuum while Table 2 in water medium. Solvation energies of these supra-molecular complexes indicate that the HCQ complexes possessed higher solvation energies with respect to their corresponding CQ complexes. The high solvation energies of the corresponding HCQ complexes revealed that in water medium, these complexes are well balanced by the hydrogen bonds. HCQ-β-CD accounts the highest binding energy value equal to −37.57 kcal/mol. The H-bonding distance that holds the HCQ moiety tightly in α–CD/β–CD hosts are in the range of 1.78–2.01 Å which clearly demonstrates that H-bonding in these supramolecular complexes are sufficiently high to offer their extra stabilities (Fig. 3 ).

Table 1

Adsorption energies; HOMO, LUMO levels, band gap and other global parameters for CQ-α-CD, CQ-β–CD, HCQ-α–CD and HCQ-β-CD inclusion complexes.

In Vacuum
	CQ	HCQ	CQ-α–CD	CQ-β-CD	HCQ-α-CD	HCQ-β-CD
E (kcal/mol)	−835,137.97	−882,501.54	−3,142,931.78	−3,527,574.99	−3,190,296.82	−3,574,942.27
ΔE (kcal/mol)			−28.22	−25.81	−29.64	−29.47
D (Debye)	5.39	6.31	9.66	10.44	15.08	12.06
HOMO (eV)	−6.97	−6.80	−7.44	−6.75	−7.33	−7.10
LUMO (eV)	−3.46	−0.10	−1.05	−0.09	−0.74	−0.42
Δ(HOMO − LUMO) (eV)	3.51	6.70	6.38	6.66	6.58	6.68
μ (eV)	5.21	3.45	4.24	3.42	4.03	3.76
χ (eV)	−5.21	−3.45	−4.24	−3.42	−4.03	−3.76
S (eV)	0.28	0.15	0.16	0.15	0.15	0.15
η (eV)	1.76	3.35	3.19	3.33	3.29	3.34
ω (eV)	7.74	1.78	2.82	1.75	2.47	2.11

Table 2

Adsorption energies; HOMO, LUMO levels band gap and other global parameters for CQ-α-CD, CQ- β–CD, HCQ-α–CD and HCQ-β-CD inclusion complexes in water medium.

In Water
	CQ	HCQ	CQ-α-CD	CQ-β-CD	HCQ-α-CD	HCQ-β-CD
E (kcal/mol)	−835,144.66	−882,511.12	−3,142,971.04	−3,527,616.81	−3,190,339.13	−3,574,985.79
ΔE (kcal/mol)			−27.44	−32.21	−31.92	−37.57
D (Debye)	7.80	8.96	11.96	16.50	20.26	17.12
HOMO (eV)	−7.00	−7.01	−7.25	−7.09	−7.06	−7.17
LUMO (eV)	−0.47	−0.44	−0.84	−0.57	−0.58	−0.63
Δ(HOMO −
LUMO) (eV)	6.54	6.57	6.41	6.52	6.47	6.54
μ (eV)	3.73	3.72	3.83	3.83	3.82	3.90
χ (eV)	−3.73	−3.72	−4.04	−3.83	−3.82	−3.90
S (eV)	0.15	0.15	0.16	0.15	0.15	0.15
η (eV)	3.27	3.28	3.21	3.26	3.24	3.27
ω (eV)	2.13	2.11	2.55	2.25	2.25	2.32

Fig. 3

Geometries of HCQ- α –CD (b) HCQ-β-CD inclusion complexes with calculated major structural parameters were listed (Angstrom unit).

Frontier molecular orbital and charge transfer characteristics

In order to demonstrate the stability of the inclusion complexes frontier molecular orbital (FMO) analysis has been done. The kinetic stability, chemical reactivity, and chemical hardness is often predicted from the energy gap between HOMO and LUMO. The global reactivity descriptors were calculated using the equations as described in the supporting information. The global reactivity descriptors are reported in Table 1 in vacuum and Table 2 in water phase. The negative chemical potential confirms that the model inclusion complexes are quite stable. The global hardness (η) of the complexes decreased after the guests penetrated into the cavity of the β-CD, while the electro negativity parameter is increasing. The calculated electrophilicity is significantly higher for the orientation A. From these results, we can conclude that Orientation A is quite electrophilic in nature. The 3D plots of the HOMO and LUMO orbitals computed at the M06-2X/6-31+G(d) level for all the complexes are illustrated in Figure S1. A careful inspection of molecular plots of encapsulated complexes, we noticed that both HOMO and LUMO contribution came from the host α–CD/β-CD and guest HCQ/CQ molecule. Such spatial ldistribution HOMO-LUMO charge densities clearly demonstrates that charge transfer involving host-guest interaction is absent in these complexes. The electrostatic potential maps (ESP) for all the complexes are plotted in the Figure S2 which confirms that for HCQ-α–CD and HCQ-β-CD complexes the dipole changes are quite prominent than their corresponding CQ complexes. We assigned stronger interactions of HCQ complexes in the α–CD/β-CD cavity is the main reason for these polarity changes. Due to this enhanced polarity changes in the HCQ, these complexes have more solvation energies than their corresponding CQ complexes.

Thermodynamic parameters analysis

From the frequency analyses different thermodynamic parameters were evaluated. The simulated thermodynamic parameters like the enthalpy change (ΔH°), the thermal Gibbs free energy (ΔG°) and entropy contribution (ΔS°) at M06–2X/6–31 G (d) levels of theory are summarized in Table 3. As illustrated from the Table 3, negative values for ΔH° and negative values for ΔG° indicated that the formation of these complexes is exothermic in nature and spontaneous at 1 atm and 298.15 K. It is quite interesting to note that all the ΔG° values are more negative for HCQ-α–CD and HCQ-β-CD than their corresponding CQ complexes. The negative enthalpy changes as the van der Waal's interaction charges form one inclusion complexes to other. The ΔS° changes for all the inclusion complexes are turned out to be negative. We note that the negative entropy change (ΔS°) is the steric barrier caused by molecular geometrical shape and the limit of α–CD and β-CD cavity to the freedom movement and rotation of guest molecule. In addition, we also note that the negative enthalpy change and entropy change for all the complexes indicated that the formation of the inclusion complex is an enthalpy-driven process and the encapsulation processes of HCQ in α–CD/β-CD host-guest complexes favorable in nature. From the adsorption, solvation energy analysis and thermodynamic parameters clearly reflects that masking of the toxic effect of CQ and HCQ along with sustainable release of drugs from the stable CQ/HCQ in α–CD/β-CD complexes is possible.

Table 3

Solvation Energy (Kcal/mol); ΔH° (Kcal); ΔG° (Kcal) and ΔS° (cal/K) for CQ-α-CD, CQ-β–CD, HCQ-α–CD and HCQ-β-CD inclusion complexes in water medium.

Orientation A	CQ-α-CD	CQ-β–CD	HCQ-α-CD	HCQ-β–CD
Solvation Energy(Kcal/mol)	−39.26	−41.82	−42.31	−43.52
ΔH°(Kcal)	−26.48	−21.96	−21.81	−24.32
ΔG°(Kcal)	−3.95	−1.16	−5.67	−4.27
ΔS°(cal/K)	−64.05	−69.984	−54.29	−67.44

Electronic absorption spectra

The UV–visible spectrum of the four supra-moleculecular complexes are studied to understand the electronic behavior of the different host guest inclusion complexes. Fig. 4 shows the simulated adsorption spectra of these inclusion complexes (data are available in Table S1). From the spectra, it is clear that all the inclusion complexes showed two major peaks at 255 nm and 326 nm. The wavelength in the vicinity of 255 nm with highest oscillator strength (0.13–0.29) corresponds to λmax which is in accordance with our previous experimental study [28]. The Table S1 indicates that there is a slight change in the effective wavelength in the vicinity of 255 nm for these complexes due to the presence of the H-bonding between the Host and the guest.

Fig. 4

Simulated UV/Vis absorption spectra for CQ-α–CD, CQ-β–CD, HCQ-α–CD, HCQ-β -CD calculated using B3LYP-D3/TZVP level of theory in water solvent.

Non covalent interactions - reduced density gradient analysis

Recently, it has been realized that NCI-RDG investigation was able to distinguish the weak interactions better than any other analysis. utilizing Non-covalent index (NCI) method we have analysed non-covalent interactions like van der Waals, H-bonding, steric communications from the plots of reduced density gradient (RDG) against electron density multiplied by the sign of the second Hessian eigen value (sign(λ2)ρ(r)) for these supra molecular complexes. Reduced density gradient is particularly helpful quantity [30] to account the non-covalent interactions. To differentiate between attractive and repulsive interactions, the eigen values (λi) of the second derivative of density is plotted within the framework of NCI technique and one can easily elucidate information about non-covalent interactions from the plots of sign(λ2)ρ vs. s. The stabilizing interactions are specified by negative values while the repulsive interaction are characterized by positive values of sign(λ2)ρ. Reduced density gradient (RDG) graphs (Fig. 5 ) indicated the presence of strong stabilizing H-bonding interactions in the HCQ/α–CD/β–CD over CQ/α–CD/β–CD inclusion complexes.

Fig. 5

Plots of reduced density gradient (RDG) against electron density multiplied by the sign of the second Hessian eigen value (sign(λ2)ρ(r)) for the (a) CQ-α–CD (b) CQ-β-CD (c) HCQ-α –CD (d) HCQ-β-CD inclusion complexes respectively.

ADMP dynamics study of the inclusion complexes

Furthermore, investigation of the stability of the inclusion complexes is of noteworthy importance toward ascertaining the validity of utilizing such macro-systems for further applications. In this context ADMP method can be considered as a reliable approach for examining the molecular stability and dynamics of macro-systems. We have examined the dynamic stabilities of the optimized geometries of for CQ-α-CD, CQ-β–CD, HCQ-α–CD and HCQ-β-CD inclusion complexes employing the ADMP method in vacuum as implemented in Gaussian09. The simulated energies are plotted within a timeframe of 300 femtosecond as illustrated in the Fig. 6 . Eventually, ADMP molecular dynamics simulations results uncover the momentous difference between the trajectories of CQ in α/β–CD and HCQ in α/β-CD inclusion complexes, where a fluctuation range of approximately of 1.0 × 10−5 and 0.2 × 10−5 Hartree for the chloroquine and hydroxychloroquine were found respectively which indicates the higher kinetic stability of the significance of HCQ in α/β-CD over CQ in α /β –CD inclusion complexes.

Fig. 6

ADMP trajectories of Etotal of 1:1 (a) CQ-α–CD, (b) HCQ-α–CD, (c) CQ-β–CD, (d) HCQ-β-CD; dynamics of optimized geometries.

The typical snapshots in Fig. 7 are chosen on three time points (0, 150, and 300 fs) confirmed that during the whole trajectory of ADMP simulations chloroquine moiety is tightly bound within the CD cavity. Fig. 8 illustrates an oscillated curve derived from the simulated trajectories of these complexes, exhibited by the variation of different O—H bonds as a function of evolution of time. Assessment of results revealed that the average H bonding distances of 1.86, 3.29, 1.96, and 1.83 Å, (for 75H-130 N, 177H-31O, 120H-128O, 114O-196H atom level as Fig. 8) have outstanding correspondence to the values in optimized geometries. Low level of fluctuations and stronger H-Bonding are the indicative of high stability of HCQ-β-CD composite systems.

Fig. 7

Typical snapshots along with the geometric parameters of for (A) CQ-α-CD, (B) CQ-β–CD, (C) HCQ-α–CD and (D) HCQ-β-CD inclusion complexesduring the ADMP dynamics simulation. Colours in graphical representations correspond to atom types as follows, red: oxygen; white: hydrogen; blue: nitrogen, and gray: carbon.

Fig. 8

.Selected H -bond distances as a function of evolution of time in the ADMP simulation process of for CQ-α-CD, CQ-β–CD, HCQ-α–CD and HCQ-β-CD inclusion complexes.

Molecular docking study of CQ and HCQ with protease of COVID-19

Furthermore, molecular docking of CQ and HCQ were done with protease of COVID-19. The binding affinities of drugs on the co-crystal of COVID-19 were evaluated by determining docking scores. We note a docking score of −6.3 for HCQ which is quite high compared to −5.9 for CQ (Fig. 9 ) which can be interpreted due to the greater number of neighboring H-bonds associated with the former than the latter.

Fig. 9

Docked structure of COVID-19 main protease with (a) CQ, (b) HCQ and their nearest neighbours (c) CQ and (d) HCQ respectively.

Conclusion

The assessment of geometric parameters, adsorption energies, and thermodynamic parameters clearly indicates that HCQ and CQ drugs formed stable host-guest inclusion complexes with α and β-Cyclodextrin. The additional H-bonding present in HCQ allows the stronger interactions of these drugs into the cavity of hosts α and β-Cyclodextrin and thereby explains the enhancement of solvation energies, adsorption energies and binding energies for these complexes. Simulation results verified that between CQ and HCQ, HCQ showed stronger affinity towards β-Cyclodextrin. Molecular docking confirmed that HCQ showed higher affinity than CQ against COVID-19 main protease. Overall the formation of these host-guest complexes is very effective strategy to mask the cytotoxicity of certain drugs and eliminating and modulating side effects of drug applications. This strategy can reduce the side effect of HCQ and CQ thereby offers a new way to use these drugs. We hope the present study should help the scientific community to develop potential therapeutics against the novel coronavirus.

CRediT authorship contribution statement

Aditi Roy: reviewed the article. Ranjoy Das: reviewed the article. Debadrita Roy: reviewed the article. Subhadeep Saha: wrote and revised the article. Narendra Nath Ghosh: did the theoretical calculations, wrote the article. Subires Bhattacharyya: reviewed the article. Mahendra Nath Roy: reviewed the article.

Declaration of Competing Interest

All the authors declare that there is no conflict of interest.