Non-competitive interactions between hydroxychloroquine and azithromycin: Systematic density functional, molecular dynamics, and docking calculations.

In this study, density functional theory (DFT) and docking calculations were systematically performed to study the non-competitive interaction between Hydroxychloroquine (HCQ) and azithromycin (AZTH). The calculated changes in Gibbs free energy and enthalpy (at 310 K) were positive, indicating the non-spontaneous formation of HCQ-AZTH specifically in water media. Docking calculation confirmed the obtained DFT result as evident from the different binding sites of both drugs to the SARS-CoV-2 main protease and human angiotensin-converting enzyme 2 (ACE2) proteins. The HCQ-AZTH structure revealed enhanced electrochemical properties, suggesting the synergy between HCQ and AZTH without affecting their therapeutic efficacy against SARS-CoV-2.

Introduction

The Corona epidemic appeared at first in Wuhan, China, in December 2019 as a new epidemic SARS-CoV-2-causing pneumonia. This corona virus disease caused by SARS-CoV-2, called COVID-19, is one of the most dangerous types of epidemics that spread and threaten human life all over the world which later become pandemic [1]. According to WHO, out of 85 million confirmed cases of COVID-19, about 1.8 million deaths were reported until now with a mortality rate of about 2.1% [2]. The danger of this disease lies in its effect on the respiratory system due to the respiratory syndrome coronavirus, which has led to a significant increase in deaths [3], [4]. So far, there is no effective treatment and preventive information related to this disease. SARS appeared about 18 years ago, and therefore the drugs that showed their effectiveness at that time may be useful to combat the current version, COVID-19 [3], [5]. Among the treatment options against COVID-19 was the trend towards reorienting old drugs such as chloroquine, hydroxychloroquine, azithromycin and others, which are well-known viral drugs with good safety profiles [3], [6], [7]. However, how Hydroxychloroquine works against COVID-19 is not yet fully clear [8], [9] and needs more practical scientific investigations and computational analysis [10], [11].

Hydroxychloroquine has long been known to be used as a treatment for malaria and other diseases. In 2020, many clinical trials were conducted to investigate the effect of hydroxychloroquine on SARS-CoV-2 patients [12]. The usefulness of using azithromycin in treating COVID-19 was also investigated [13]. One of the most prominent research results was the observation of a positive effect using hydroxychloroquine in combination (one at a time or both together) with azithromycin [9], [14]. Thus, it is recommended to use a combined therapeutic dose of hydroxychloroquine and azithromycin to treat COVID-19 and prevent its spread. However, the association of these protocols with deaths from thrombosis is not entirely clear, especially in patients who have undergone a prolonged course of treatment. Therefore, the research related to these two drugs is still ongoing and of great importance, including what we will discuss in our current study regarding the interaction between the two drugs and their synergy in treating COVID-19. The possibility of HCQ-AZTH interaction arising from their co-administration for COVID-19 treatment protocol cannot be ruled out as pointed out by Gaisenok [15]. Studying this interaction is therefore very important for both communities of researchers who are supporting or discouraging the HCQ-AZTH treatment protocols.

Azithromycin is one of the most popular effective antibiotics. Among the reasons that make azithromycin suitable as a drug within the protocols for treating Corona is its long-standing use as a treatment for respiratory diseases and its ability to reduce the spread of bacterial infection. This is in addition to the scientific reports that showed its ability to interact with viral enzymes, such as those produced by the Corona virus [13]. The SARS-CoV-2 main protease 6LU7 (Hereafter to be referred to as Mpro) is a key enzyme used as a therapeutic target for preventing viral replication [16], [17]. Also, the angiotensin-converting enzyme 2 (ACE2) is the main human cell receptor required for SARS-CoV-2 viral entry [18]. Therefore, we assessed the antagonistic effects of the individual use of HCQ and AZTH as well as combined HCQ-AZTH by their binding and targeting Mpro and ACE2 proteins to fight SARS-CoV-2.

From the above, the effectiveness of the two drugs (hydroxychloroquine and azithromycin) is shown separately in combating COVID-19. It remains necessary to study the effect of the interaction between the two drugs on their therapeutic properties, especially against the SARS-CoV-2.This study will also contribute to the renowned and active Pharmacokinetics field of research concerned with drug-drug interaction [19].We have conducted this study using computational methods and found that the two drugs cannot interact spontaneously at the human body temperature, despite their ability to form a thermodynamically unstable composition Thus, our results show that the interaction between these two drugs does not change the therapeutic properties of each separately. We have verified this by performing virtual screening (docking) calculations to study the locations of the drugs' targeting human and viral proteins specific to SARS-CoV-2. Significantly, we found that both drugs targeted different sites on these proteins, working as non-competitive antagonists in treating COVID-19 disease. Our study is pertinent to the previous study on the role of interaction between the two drugs in compacting resistant malaria [20]. This study determined that there is no evidence that the interactions between the two drugs are related to malaria-resistance treatment. Moreover, and as an extension to the present study, we will consider in our future investigations the possibility of AZTH-HCQ co-crystal [21] formation which improves the drugs’ physicochemical properties while maintaining their pharmacological characteristics.

In this work, we demonstrated two important points: i) HCQ and AZTH do not seem to form stable complexes in water solution at a temperature close to the human body (about 37 °C) and ii) we confirmed that HCQ and AZTH have different binding sites on both the SARS-CoV-2 main protease and the human angiotensin-converting enzyme, the two possible targets for the COVID-19 treatment. Additionally, our calculation showed enhancements in the AZTH-HCQ structure in terms of its higher polarizability, the smaller energy gap between its frontiers orbitals, and other electrochemical properties as compared to those properties for the individual AZTH and HCQ structures [10].

Materials and methods

The crystal structures of SARS-CoV-2 main protease, Mpro, and human cell receptor for SARS-CoV-2, ACE2, were obtained from the protein-databank (PDB)[22] with PDB codes 6LU7 and 1R42, respectively.

The structures of the drug ligands, hydroxychloroquine (C18H26ClN3O) and azithromycin (C38H72N2O12), were obtained from PubChem (pubchem.ncbi.nlm.nih.gov/) with PubChem CIDs as 3652, and 447043, respectively. Interconversion of chemical structures between different formats, as required by specific programs in this study, was performed by Open Babel software [23].

The calculations concerned with the polarizability and the optimization of the different conformation of HCQ-AZTH in the gas phase and water media have been carried out using the B3LYP functional and def2-TZVP basis set with D3BJ dispersion correction, as implemented in the ORCA program [24]. The hybrid B3LYP functional is known for its accuracy and is the most popular functional in chemistry as pointed out by Zhang et al [25]. It is usually used for calculating the electronic properties of organic compounds (e.g. [26]). Our results (e.g. the polarizability) are thus consistent with existing reports (e.g. [10]). While decent, our def2-TZVP basis set is not large enough to eliminate basis set superposition error (BSSE), which will lead to spurious overbending. Therefore, a counterpoise calculation was attempted to ascertain the amount of BSSE. The counterpoise correction to the binding energy was performed following the procedure described in the ORCA manual (version 4.2.1) Section 8.1.6. Briefly, the main part of this procedure is concerned with equation 8.10 in this manual which gives the BSSE from the energies of AZTH-HCQ, HCQ and AZTH calculated at their respective optimized geometries and the optimized geometry of AZTH-HCQ. For the vibrational frequency calculation (Infrared spectra) we use the numerical frequency mode, which is computationally expensive, yet informative, compared to the analytical frequency one. The calculations in the water media were obtained by the conductor polarizable continuum model (CPCM) [27].The molecular electrostatic potential (ESP) was calculated using IQmol software [28].

Regarding the calculations of the vibrational frequencies, we observed no negative(imaginary) frequencies, suggesting the mechanical stability of our optimized structures. For all structures the charge was zero and multiplicity (M = 2S + 1, where S is the total spin) was 1.

Protein-Ligand interaction was systematically analyzed by the AutoDock-vina docking program [29]. Water molecules were removed from the target proteins using UCSF Chimera software [30] in preparation for the subsequent docking stage. We employ the setting where the ligands (compounds) were flexible and the receptors (protease) were rigid. The size of search space for the best binding sites was maximized (~50 Å ×65 Å× 65 Å search space grid with ~ 1 Å spacing) to include the whole protein surface. The number of docking configurations was 10. Prior to the docking process, the structures of the ligands underwent energy minimization using the steepest-decent-algorithm as implemented in PyRx (Python Prescription Virtual Screening) module [31]. Several orientations (~10–12 orientations) of the ligands were searched, and the configuration with the best docking score was selected. The protein–ligand 2-dimensional interaction plots were generated using the Biovia Discovery Studio [32].Visualizations of the optimized and docked structures were performed with IQmol, Avogadro [33] or Discovery Studio.

The molecular dynamic (MD) [34] calculations were performed on the best-docked protein–ligand complex structures using NAMD (abbreviation for Nanoscale Molecular Dynamics) software [35] and VMD (Visual Molecular Dynamics package [36]) to visualize and analyze the molecular dynamic results. The best-docked ligand structure in “pdbqt” AutoDock output format was converted into “pdb” format using Openbabel [37] utility and then properly combined with the protein (Mpro or 1R42 (ACE2)) “pdb” file to form the protein–ligand complex structure needed for the MD simulations. VMD was then used to extract protein and ligand coordinates from the complex structure file and to generate the protein-structure-file “psf” for the protein. The simulation utilized the CHARMM (Chemistry at Harvard Macromolecular Mechanics) General Force Field (CGenFF) [38], [39] that can be imposed on the ligand from the CHARMM-GUI web service [40]. This produced the protein-structure-file “psf” with information about the force terms and ligand geometry required by NAMD. Before the CGenFF step, the ligand structure was protonated (adding hydrogen) and converted from “pdb” to “mol2” format (suitable for CGenFF procedure) using Openbabel module.

The “psf” and “pdb” files thus produced for the protein and ligand were used by a suitable VMD script to generate a combined complex structure in “psf” and “pdb” formats that were needed for the solvation step by VMD. The complex was solvated with a water box using the complex’s structure dimension.

We begin by specifying the periodic boundary condition for the solvated protein–ligand complex using VMD. The execution of the MD calculation utilized 1000 steps of structure optimization followed by 200,000 MD steps with 1 fs (1 × 10−15 s) per step (i.e 0.2 ns MD calculations). The calculation temperature was 300 K. The output “par” files (from CHARMM GUI step), containing information about the topology and force field for water (using TIP3P [transferable intermolecular potential with 3 points] water potential [41]) and other chemical materials, were provided for the NAMD calculations. Default NAMD settings were used for the remaining parameters. We restrict ourselves to the analysis of the RMSD (Root-Mean-Square-Deviation) for the complex conformation using VMD.

Results

Molecular electrostatic potential (ESP) and infrared spectra calculations

The electrostatic potentials (ESP) isosurfaces of individual HCQ and AZTH are, respectively, shown in Fig. 1 (a) and 1(b). The red and blue isosurfaces correspond to more negative and positive ESP, respectively. The color bar shows the variation of ESP between these limits. Fig. 2 (a) is the structure for the interacting AZTH-HCQ initially formed manually by bringing the most negative HCQ potential in proximity with the most positive AZTH potential. Then the structure was geometrically optimized giving rise to the binding energy of – 18.68 kcal/mol (–0.81 eV). Another AZTH-HCQ configuration with the binding energy of 376.35 kcal/mol (~16.32 eV) is shown in Fig. 2(b). Due to symmetry, other configurations are expected to give rise to binding energy greater than that in Fig. 2(a). Thus, the structure of Fig. 2(a) has been considered for all the subsequent investigations in this work. We followed a similar procedure as that reported by Novir et al [42] for selecting the AZTH-HCQ orientation for optimum interaction energy.

Fig. 1

Calculations in a gas phase of the electrostatic potential isosurfaces (EPS) for the optimized AZTH (a) and HCQ (b) structures where the red and blue represent most negative and most positive EPS.

Fig. 2

(a) The configuration for the structurally optimized AZTH-HCQ complex interacting through their most negative EPS and most positive EPS sites as shown in Fig. 1. The binding energy for this configuration is −18.68 kcal/mol (b) another AZTH-HCQ configuration with binding energy 376.35 kcal/mol.

Briefly, we tackled the computationally demanding geometry optimization of different AZTH-HCQ structures (each with > 180 atoms) in the quest for the one with the least binding energy. The structures span all possible orientation of HCQ around AZTH. The selected structure has a charge distribution corresponding to the electrostatic-potential profile shown in Fig. 1.

Frequency calculation was conducted to verify the stability of the optimized HCQ, AZTH, and AZTH-HCQ (shown in Fig. 2(a)). Fig. 3 shows the infrared spectra for these structures in water media (left panels) and gas phase (right panels). Clearly, no imaginary (negative) frequency modes were obtained for these spectra, strongly suggesting that the local minima optimizations were achieved. The correction in the vibrational entropy was incorporated in the calculation of the Gibbs free energy following Grimme et al [43] method of the dispersion-corrected density functional as implemented in ORCA software. The entropy corrections in gas phase and water media calculations were –5.059 eV and –5.056 eV, respectively, with the latter slightly less in magnitude than the former. A more rigorous analysis of the entropy as reported by Falivene et al [44] may lead to improved values of the entropy corrections for better comparison with experimental results. However, our present analysis may be sufficient at least for qualitative discussions. Specifically, our data for the Gibbs free energies for AZTH and HCQ are consistent with the reported ones (e.g. [10]).

Fig. 3

Calculated Infrared spectra (calculated with the numerical frequency mode) for HCQ, AZTH, and AZTH-HCQ in both gas phase and water medium.

Electrochemical properties of the structures

The frontier orbitals (Highest Occupied Molecular Orbitals (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO)) are important for assessing the drug properties. It has been reported that the degree of the reactivity of the drug can be associated with the energy gap between HOMO and LUMO [45], [46]. From these descriptors, we can get several important electrochemical parameters such as the chemical hardness, η, electrophilicity, ω, etc … as listed for the structures (AZTH, HCQ, and AZTH-HCQ) in Table 1 for the calculations in both gas phase and water medium. The energy gap is defined as E g= (ELUMO-EHOMO). Fig. 6 shows the isosurfaces for the HOMO and LUMO orbitals for the AZTH-HCQ structure in water media. No significant changes were visualized for the isosurfaces in the gas phase (data not shown). The calculated polarizabilities, α, and dipole moments, μ, for all the structures are also listed in Table 1 in the convenient units of Bohr3, and Debye, respectively.

Table 1

Calculated thermoelectronic energies (at the normal human body temperature of 310 K) for AZTH, HCQ, and AZTH-HCQ in gas phase and water media.

	Gas			Water
	AZTH	HCQ	AZTH-HCQ	AZTH	HCQ	AZTH-HCQ
Eb (kcal/mol)	–	–	– 18.68	–	–	– 12.00
ΔG (eV)	–	–	– 0.15	–	–	+0.15
ΔH (eV)	–	–	– 0.18	–	–	+0.12
EHOMO (eV)	−4.95	−5.54	−4.79	−5.18	−5.58	−5.23
ELUMO (eV)	0.01	−1.14	−1.27	−0.12	−1.30	−1.29
Eg (eV)	4.96	4.41	3.52	5.06	4.28	3.95
χ (eV)	2.47	3.34	3.03	2.65	3.44	3.26
η (eV)	2.48	2.20	1.76	2.53	2.14	1.97
ω (eV)	1.23	2.53	2.60	1.39	2.77	2.69
ECT	–	–	−0.52	–	–	−0.56
μ (Debye)	4.83	6.08	7.04	7.23	8.77	10.17
α (Bohr3)	467.35	235.43	704.45	638.56	335.57	983.33

Fig. 6

The calculated HOMO and LUMO isosurfaces for the AZTH-HCQ complex structure in water phase. The value of the plotted isosurfaces is 0.02000 a.u.

Docking AZTH and HCQ on SARS-CoV and human receptors

To gain more insight into the mechanism of interaction between the drugs (HCQ and AZTH) and SARS-CoV-2main protease (Mpro) as well as the human angiotensin-converting enzyme 2 (ACE2 or 1R42), molecular docking was conducted using Autodock vina software. The docking results are shown in Fig. 4 (a) and 4(b) which illustrate the 3-dimensional (left panel) and 2-dimensional (middle and right panels) binding sites of the drug ligands (AZTH and HCQ) on Mpro and ACE2 receptors, respectively. The abbreviations within the colored circles (2-dimensional interaction of Fig. 4) represent the amino acid residues for Mpro and ACE2 involved in the interaction. For instance, the Mpro residues, ASN:238 and LYS:137, form two hydrogen bonds to the AZTH (dark green dotted lines in the middle panel of Fig. 4(a)). Similarly, three hydrogen bonds are formed between the AZTH and the ACE2 residues ASP:350, PHE:390, and ARG:514 (dark green dotted lines in the middle panel of Fig. 4(b)). The calculated docking affinities (or energies) are −7.0, −9.4, −5.7, and −6.1 kcal/mol for AZTH-Mpro, AZTH-ACE2, HCQ-Mpro, and HCQ-ACE2 interactions, respectively, as shown in Fig. 4.

Fig. 4

Calculated docking interactions of AZTH and HCQ with (a) Mpro protease and (b) ACE2 (or 1R42) protein. The left colored sheets and helixes are 3-dimensional models for the proteins-drugs interactions. The docking energies (affinities) are shown in kcal/mol beneath the corresponding 2D-interaction models shown to the right of the 3-dimensional model. The colored circles in the 2D models show the names of the protein amino acid residues within the interaction pockets.

Molecular dynamic (MD) calculations

Fig. 5 illustrates the MD Root-Mean-Square-Deviation (RMSD) versus the time (in nano-s, ns) for the atomic movements in the receptor-ligand complex systems corresponding to the best-docked structures shown in Fig. 4 where Mpro and 1R42 (ACE2) represent the receptors, and HCQ and AZTH are the ligands.

Fig. 5

The protein–ligand Root-Mean-Square-Deviation (RMSD) as calculated during 0.2 × 10−9 s (0.2 ns) molecular dynamic simulations. The proteins are 6LU7 (or Mpro) and 1R42 (or ACE2), and the ligands are AZTH and HCQ.

Discussion

As we mentioned above, the AZTH-HCQ structure of Fig. 2(a) was selected based on the ESP profile (Fig. 1) and its lowest binding energy, E b ~ – 18.68 kcal/mol (–0.81 eV) in the gas phase (Table 1). The calculated infrared spectra (Fig. 3) for this structure and its individual components (AZTH and HCQ) exhibit no negative frequencies, suggesting the achievement of the local minima (structural stability) in our geometry optimizations. For instance, the optimized C-C bond length for HCQ is ~ 1.418 Å (gas phase), agreeing with the theoretical and experimental reports [47]. This C—C bond length has slightly reduced to 1.416 Å for AZTH-HCQ. We emphasize that by C—C bond length we mean a length between two nearest C neighbors in one of the HCQ rings structures, that was selected for a rapid assessment of the optimized geometrical structure.

In Fig. 3, the transmission valleys (or absorption peaks) around frequencies (or wavenumbers) > ~ 3000 cm−1 is associated with O—H, N—H, and C—H bonds available for most of the drugs. While for < 2000 cm−1 the transmission valleys are mostly due to absorption of the spectrum by CC and CN type of bonds. Comparing the spectra in the gas phase and water media reveals that the effect of the water is to enhance the transmittance as compared to that in the gas phase. In general, no new peak has appeared in the AZTH-HCQ spectra when compared to that for the individual HCQ and AZTH compounds, which might indicate the weak AZTH-HCQ coupling as will be supported by the DFT and docking calculations.

The consequence of this type of coupling (or interaction) can be made evident from the calculated (at 310 K) changes in Gibbs free energy (ΔG) and enthalpy (ΔH) shown in Table 1. Clearly, both changes are positive, with those for the gas phase being lower than the water medium. The positive ΔG and ΔH suggest no spontaneous formation of AZTH-HCQ complex specifically in the water phase, yet, the complex is expected to exist at temperatures > 310 K at thermodynamically unstable states. Additionally, the calculated binding energy from DFT is E b = E complex − (E HCQ + E AZTH) − E BSSE, where E complex, E AZTH and E HCQ are the calculated energies for AZTH-HCQ complex, AZTH, and HCQ, respectively. The Basis set superposition error (BSSE) was estimated as – 0.18 eV from counterpoise calculation [48] in water medium. The same correction value of – 0.18 eV was used for the gas phase structure because the optimized atomic sites in the two calculations relating to the gas phase and water medium are almost similar. From Table 1, the binding energy in the water medium is attenuated in comparison to that in the gas phase. This may be associated with the dielectric properties of the water media. Therefore, AZTH-HCQ interaction does not significantly alter the pharmacokinetics and therapeutic properties of the individual components (AZTH and HCQ). This may act as a support for the existing experimental report [20] which excludes that the interaction between the two drugs is the reason for the success of the treatment of chloroquine-resistant malaria.

To further emphasize the non-competitive nature of the interaction between AZTH and HCQ, we have calculated the change in Gibbs free energy, ΔG, and enthalpy, ΔH, at the room temperature of 298 K. Our results show that the interaction is thermodynamically favorable with ΔH < 0 and ΔG < 0 for both gas-phase and water medium cases. This indicates that the interaction between the two drugs is attenuated within the human body. It may also suggest that the two drugs must not be mixed in a solution at room temperature for oral administration to patients with swallowing difficulties.

The above-mentioned nature of the interaction between AZTH and HCQ stimulated us to study (via molecular docking [49] calculations) the ability of the two drugs to target the Corona virus and other human proteins in a non-competitive manner.

The molecular docking gives us more insight into the binding mechanism of the drug ligands (HCQ and AZTH) on the SARS-CoV-2 protease and human protein receptors. The docking results of HCQ-Mpro (Fig. 4(a)) indicate a binding affinity of −5.7 kcal/mol. Our calculated binding energy and binding site (e.g. around the residues HIS-41, MET-49, MET-165, GLU-166, LEU-167, ASP-187, and GLN-192 for HCQ-Mpro in Fig. 4(a)) agree very well with those obtained by Hagar et al [26]. The main essence of Fig. 4 is to show the non-competitive screening capabilities of HCQ and AZTH against Sars-Cov2 main protease (Mpro) and the human angiotensin-converting enzyme 2 (ACE2 or 1R42) as evident from the different interaction pockets for each drug on these receptors (Mpro and ACE2). The left 3-dimensional panels in Fig. 4 clearly indicate the different interaction pockets for the drugs (structures dominated by grey bars) within Mpro and ACE2 (shown as colored sheets and helixes). The corresponding amino acid residues involved in the ligand–protein interaction are depicted by the 2-dimensional models in the middle and right panels of Fig. 4. Here, HCQ and AZTH strongly bind to the ACE2 (Fig. 4(b)) with −6.1 and −9.4 kcal/mol binding affinities, respectively, as compared with their binding to the Mpro (Fig. 4(a)). One of the merits of the ligand–protein binding (either on ACE2 or 1R42) is to hinder the viral replication process. For example, the drug-ACE2 binding would modify the morphology of ACE2, reducing the chance of SARS-CoV-2 replication process which utilizes ACE2 as a cell coupling receptor [50].

To support the docking results, we have performed a systematic Molecular-Dynamic calculation (Fig. 5) on the best-docked (see Fig. 4) complex structures. The RMSDs exhibit nearly plateaus regions which may be regarded as manifestations for the equilibrium states reached within around 1 ps (1 pico-s). The average RMSD for these systems is ~ 2 Å, with HCQ-Mpro and AZTH-Mpro complexes attaining the least RMSD. As implied for example by Beura et al [51] RMSD < 3 Å is an acceptable indicator for the conformational stability of small ligand–protein complexes. We expect the equilibrium states (plateau regions) in Fig. 5 to be preserved for longer than 0.2 ns. The effect of the solvation layer (e.g. water in our case) on the stability of the shell and core protein structure during the MD calculation is an important factor for MD simulation. It is mainly pertinent to the stability of macromolecules which requires an extended time ~ nanoseconds of MD calculations to study the equilibration of the system. Nevertheless, our present simulation time of 0.2 ns may be enough to validate the stability of the docking results of the small ligands in the best-docked ligand-receptor (AZTH-protein or HCQ-protein, Fig. 4) configurations.

We now discuss the calculated electrochemical properties (in gas and water conditions as shown in Table 1) for the drug. Clearly, the energy gap, E g = E LUMO– E HOMO, between the HOMO and LUMO is reduced for the AZTH-HCQ as compared to that for the isolated AZTH and HCQ. Smaller E g means softer structures with smaller chemical hardness, η = –(E HOMO − E LUMO)/2, as can be seen from Table 1. Drugs with smaller η are pharmaceutically favorable due to their ability to exchange charges with the target proteins [45], [46], [52], [53]. We can involve η in some quantitative expressions related to the electrophilicity ω = μ2/2η, where the chemical potential μ = (E HOMO + E LUMO)/2 from which the electronegativity is defined as χ = – μ. An important descriptor is the electrophilicity-based charge transfer defined as ECT = (ω/2χ)AZTH – (ω/2χ)HCQ, such that if ECT > 0, AZTH acts as an electron acceptor and if ECT < 0, AZTH acts as an electron donor. From ECT values in Table 1 we see that AZTH donates electrons to HCQ. This may nicely agree with the HOMO and LUMO distributions over the AZTH and HCQ, respectively, seen in Fig. 6 .

Moreover, the calculated polarizabilities (α) and dipole moments (μ) are significantly enhanced for the complex structures as compared with the noninteracting ones (Table 1). As expected, the more pronounced values are obtained in water media. We also notice that the complex polarizabilities are nearly the additive of those for the noninteracting compounds, indicating the weaker AZTH-HCQ coupling. Our data for α and μ well conform with the reported ones for individual HCQ and AZTH [10].

Conclusion

We have managed to assess the pharmacokinetics modification of the azithromycin (AZTH) and hydroxychloroquine (HCQ) that could have emerged due to drug-drug interactions. This was accomplished through systematic and consistent DFT and virtual screening (docking) calculations. For the combined AZTH-HCQ we achieved positive changes in the Gibbs free energy and enthalpy (at human body temperature of 310 K), indicating a nonspontaneous drug-drug interaction and highlighting the safe combinational administration of both drugs. Noncompetitive targeting of viral or human proteins by these drugs is an important consequence of their safe interaction which was verified through molecular docking calculations. Our results validate the synergy of HCQ and AZTH when given in combination without affecting their individual therapeutic efficacy against SARS-CoV-2, agreeing with the existing experimental report. In other words, we demonstrated two important points: i) HCQ and AZTH do not seem to form stable complexes in water solution at a temperature close to the human body (about 37 °C) and ii) we confirmed that HCQ and AZTH have different binding sites on both the SARS-CoV-2 main protease and the human angiotensin-converting enzyme, the two possible targets for the COVID-19 treatment.

Funding

This work was supported by The Fast Track Funding Path for Coronavirus (COVID-19) project funded by The King Abdulaziz City for Science and Technology (KACST) underGrant No 5–20-01–006-0021

CRediT authorship contribution statement

Mohammed A.H. Khalafalla: Conceptualization, Formal analysis, Methodology, Project administration, Software, Visualization, Writing - original draft, Writing - review & editing. Chokri Hadj Belgacem: Investigation, Methodology, Resources, Writing - review & editing. Ismail Abdelrehim: Investigation, Validation, Writing - review & editing. Kamel Chaieb: Investigation, Methodology, Writing - original draft.

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.