DFT and molecular docking study of chloroquine derivatives as antiviral to coronavirus COVID-19

2.1 DFT calculations

The GaussView program (GaussView, Guassian, Inc.) was utilized to model the initial structures of the chloroquine and the chloroquine phosphate molecules. Subsequently, their molecular geometries optimizations were carried out in the gas phase with the density functional theory (DFT) with the Gaussian 09 software package (Gaussian 09, Revision C.01, Frisch et al., 2009). All the quantum-chemical calculations have been performed via the hybrid B3LYP (Becke’s three parameter hybrid functional with Lee-Yang-Parr correlation functional LYP (Lee et al., 1988; Becke, 1993) at 6-31G* basis set. Furthermore, several electronic properties for instance the frontier molecular orbitals, gap energies, reactivity descriptors were computed using TD-DFT approach (Liu et al., 2015; Becke, 1993). The density of states (DOS) plots was obtained by using Gauss-Sum software (O'Boyle et al., 2008).

2.2 Ligands and proteins preparation

The 3D structures of COVID-19 protein were retrieved from the RCSB PDB database (http://www.rcsb.org) (http://www.rcsb.org/). The Protein Data Bank (PDB) archive contains thousand protein structures obtained either by crystallography X-ray or by NMR. Concerning ligands, the 2D structures of chloroquine and chloroquine phosphate were extracted from the PubChem online database (https://pubchem.ncbi.nlm.nih.gov/). The ligands were saved in the MDL Mol file format. Then, they were converted to a PDB file format by using Accelrys Discovery Studio Visualizer (Visualizer, 2005). Thereafter, Rapid-Screening docking was carried out using iGEMDOCK program (Yang and Chen, 2004). It is a Drug Design System for docking calculations and screening by BioXGEM labs. All the trials were docked with a population size set to 800, with 80 generations and 10 solutions.

3.1 Optimization of chloroquine and chloroquine phosphate

Optimized structures and numbering of atoms of chloroquine and chloroquine phosphate molecules are shown graphically in Figs. 1 and 2, obtained at B3LYP/6-31G* method. Table 1 illustrates their geometrical parameters such as the calculated total energies, the dipole moments, the RMS and the maximum Cartesian force. The global minimum energies are found to be −1326.0352 a.u (≈ −36083 eV) and −2614.3242 a.u (≈ −71139) for chloroquine and chloroquine phosphate, respectively. The RMS Cartesian force values are equal to 2.412 0.10⁻⁶, 0.04067 in chloroquine and chloroquine phosphate. Their maximum Cartesian forces are found to be 8.593 0.10⁻⁶ and 0.1449. The dipole moment of a molecule is given in the form of a three-dimensional vector and which reflects the molecular charge distribution. Hence, it can be employed as a descriptor to describe the charge movement throughout the molecule. As a result of DFT/B3LYP/6-31G* calculations, the highest dipole moment was observed for the chloroquine phosphate (∼24.49 Debye) whereas the smallest one was observed for the chloroquine (∼6.05 Debye). Of course, the adding of other atoms in the geometry of the chloroquine has an influence on their stability. We can notice that the chloroquine compound becomes more stable when adding the phosphate groups since the global minimum energy decreases. Also, the strong increase in the dipole moment value shows that the chloroquine is harder before adding the phosphate groups. Moreover, it promotes the formation of hydrogen bonds.

Close

Fig. 1

Optimized structure of the chloroquine by using DFT/B3LYP/6-31G* method.

Export to PPT

Close

Fig. 2

Optimized structure of the chloroquine phosphate molecule.

Export to PPT

The optimized geometrical parameters of chloroquine derivatives have been determined by the above method and they are given in Tables 2 and 3 with the experimental bond angles and bond lengths. First, we observed that the theoretical bond lengths of chloroquine compound are almost similar with the experimental results (Busetta and Courseille, 1973), since the value of RMSD is very small (0.001 Å). The same applies to the bond angles which have an RMSD value equal to 0.298°. Same thing for the chloroquine phosphate, according to the result as collected in table 3 the bond distances and bond angles show good agreement with the experimental data (Albesa-Jové et al., 2008). We find that the RMSD value is equal to 0.065 Å for the bond distances and 3.382° for the bond angles. Results reveal that the carbon–carbon bond distances are found in the range 1.374–1.546 Å for C₂₀-C₂₂ and C₅-C₇, respectively for the chloroquine. In the benzene ring (I), the carbon–carbon bond lengths C₁₃-C₁₇, C₁₃-C₁₈, C₁₇-C₂₀, C₁₈-C₂₁, C₂₀-C₂₂ and C₂₁-C₂₂ are 1.435, 1.418, 1.421, 1.378, 1.374 and 1.411 Å, respectively. The C–C bond alienation in the pyridine ring (II) is between 1.394 Å (for C₁₂-C₁₆ bond) and 1.445 Å (for C₁₂-C₁₃ bond). While, for chloroquine phosphate, the bond length between two carbon–carbon in the two rings is in the range 1.383–1.419 Å for benzene and 1.366 to 1.464 Å for pyridine ring. It is seen that the B3LYP calculated hydrogen bonding distances C–H vary from 1.009 Å (for N₃-H₃₀) to 1.099 Å (for C₅-H₂₄) for chloroquine and from 1.084 Å (for C₁₀-H₂₇ bond) to 1.524 Å (for C₂₁-C₂₂ bond) for chloroquine phosphate. Three nitrogen N atoms exist in the structure of chloroquine: the order of the N-C bond length is N₂-C₁₀ > N₂-C₁₁ > N₂-C₈ > N₃-C₇ > N₃-C₁₂ > N₄-C₁₇ > N₄-C₁₉ having values 1.470 > 1.469 > 1.467 > 1.465 > 1.370 > 1.365 > 1.319 Å, respectively. The bond distance of N₃-H₃₀ is equal to 1.009 Å. The bond angle of chloroquine between the C₇-N₃-H₃₀ and C₁₂-N₃-H₃₀ are ∼ 115.047° and ∼ 116.505°, respectively. Concerning the chloroquine phosphate, we note that the single N₅-C₆ bond length of 1.387 Å for ring pyridine is higher than the N₅-C₄ double bond (1.353 Å). The P-O bond lengths are obtained to be in range 1.489–1.693 Å (for P₅₈-O₆₁ and P₅₈-O₆₂). The O-P-O bond angles are reported in range 107.7–112.02°, whereas it is computed in range 102.543–124.278°. The C₈-Cl bond length is observed at 1.743 Å and calculated at 1.748 Å. The C₉-C₈-Cl and C₈-C₉-C₁₀ bond angles are at 119.733° and 116.940°, respectively.

3.2 Frontier orbitals and quantum chemical calculations

Frontier molecular orbitals (FMOs) often play dominant roles in molecular systems. The fundamental idea of this theory can be abridged in the form of a simple rule telling the condition for a simple course of the reaction by the requirement of the maximal positive overlap between LUMO (empty state) and HOMO (filled state) orbitals. LUMO (lowest unoccupied molecular orbital) is directly related to electron affinity, while HOMO (highest occupied molecular orbital) is related to ionization potential (Xavier and Periandy, 2015; Abraham et al., 2017). These orbitals help to understand the chemical stability and the reactivity of the molecule (Asiri et al., 2011; Kosar, 2011). In order to predict the energetic behaviors and the reactivity of the chloroquine and the chloroquine phosphate against COVID-19 virus, the FMOs in the electronic transitions and their energies difference Eg are determined. A detailed analysis of the HOMOs and LUMOs orbitals is listed in Table 4, where orbital energies, energy band gap and reactivity descriptors (like electron affinity, chemical softness, ionization potential, chemical softness….) are reported. The gap between two energetic states describes the molecular chemical reactivity. The energies of the four important FMOs (HOMO, HOMO − 1, LUMO and LUMO + 1) were calculated via the TD-DFT approach with B3LYP/6-31G* level. Their 3D plots are illustrated in Figs. 3 and 4. It is clear from the figure of the chloroquine molecule that the HOMO and LUMO orbitals are localized essentially on the benzene and pyridine rings. The green color represents the negative phase; on the other hand the red color corresponds to the positive phase which is well clarified in the density of states (DOS) spectrum (Fig. 5). DOS spectrums characterize the energy levels per unit energy increment and its composing in energy. The displaying study per orbital shows that the green and the red lines in these curves correspond to the HOMO and LUMO energy levels, respectively. As a result, the energy level of the HOMO orbital is about −5.594 eV and the energy level of the LUMO orbital is about −1.115 eV. The HOMO-LUMO gap energy (Eg) of the chloroquine is equal to −4.479 eV. This low energy value promotes the transfer of electrons in the chloroquine molecule. These values are compatible with those obtained by the DOS spectrum. The state HOMO-1 form another set of degenerate orbital −5.747 eV lower in energy than the HOMO set. As shown for the chloroquine phosphate, LUMO orbital lying at −2.59 eV, located on all the atoms of the benzene and pyridine rings. The HOMO orbital is lying at −5.228 eV. Consequently, Eg is closed to −2.629 eV. The change observed here in the gap value from −4.479 eV to −2.629 eV in solution involves an expected high reactivity for the chloroquine phosphate. This decrease in gap energy makes the flow of electrons easier, so the molecule becomes soft and more reactive. We can also note that the chloroquine molecule is harder before adding the phosphate groups, given the energy value of gap. This result is in agreement with the strong increase in the dipole moment value of 6.05 Debye (of chloroquine) to 24.49 Debye (of chloroquine phosphate).

Close

Fig. 3

The atomic orbital compositions of the HOMO, HOMO-1, LUMO and LUMO + 1 frontier molecular orbitals for chloroquine molecule.

Export to PPT

Close

Fig. 4

The atomic orbital compositions of the HOMO, HOMO-1, LUMO and LUMO + 1 frontier molecular orbitals for chloroquine phosphate.

Export to PPT

Close

Fig. 5

DOS spectrum of chloroquine (a) and chloroquine phosphate (b) molecules.

Export to PPT

Using the energies of FMOs, we calculated the reactivity descriptors of chloroquine and chloroquine phosphate molecules. A =  − E_LUMO: represent the electron affinity; I =  − E_HOMO represent the ionization potential and μ = 1/2(I + A) is the electronic chemical potential. The chemical hardness (η) is found to be 2.239 and 2.629 eV for chloroquine and chloroquine phosphate, respectively. The chemical softness (ζ) has been computed and found to be 1.1195 and 1.3145 eV⁻¹. Moreover, the electrophilicity index (ω) is about 2.512 eV for chlroquine and 2.912 eV for chloroquine phosphate. Based on the value found of the electrophilicity index, we can conclude that the chloroquine phosphate is a good electrophile better than chloroquine. Therefore, it is able to accept an electron doublet in order to form bonds with another reagent which is necessarily a nucleophile. Electronegativity is also determined (χ = (I + A)/2) and it is found to be χ _chloroquine = 3.3545 eV and χ _{chloroquine phosphate} = 3.9135 eV.

3.3 Molecular electrostatic potential

The molecular electrostatic potential (MEP) is a well-established tool for the study of molecular reactive properties and to describe intermolecular interactions (Reed and Weinhold, 1985). It allows us to search the most reactive nucleophilic and electrophilic sites of a molecule against the reactive biological potentials (Gökce et al., 2013). These sites promote the formation of hydrogen bonds. The electrophilic site indicates a strong attraction, while the nucleophilic site indicates a strong repulsion. The electrostatic potential diagrams of chloroquine and chloroquine phosphate are illustrated in Fig. 6 at B3LYP/6-31G* method. MEP diagram gives negative, positive and neutral electrostatic potential regions in terms of color grading and is an indicator in the research of molecular structure properties. The red color represents the most electronegative electrostatic potential. That is, atoms in this region have a tendency to attract electrons (electrophilic). The blue color indicates the most electropositive potential (strong attraction) and the red color indicates the most electronegative potential (strong repulsion). Regions where the potentials are zero are denoted by green color. As a results, MEP surfaces varies between −5.504 0.10⁻² a.u (deepest red) to 5.504 0.10⁻² a.u (deepest blue) for chloroquine and between −0.116 a.u to 0.116 a.u for chloroquine phosphate. As can be seen, the MEP map of chloroquine molecule (Fig. 6a), a maximum positive region is localized on the nitrogen N₃ and hydrogen H₃₀ atoms indicating a possible site for electrophilic attack. The zero potential sites (green color) are found in the benzene ring. For the chloroquine phosphate (Fig. 6b), the positive potential (blue and light blue) sites are found in the benzene and pyridine rings (electrophilic reactivity). It can be inferred that the oxygen atoms O₆₁ and O₆₂ indicate the neutral potential of the molecule.

Close

Fig. 6

Molecular electrostatic potential (MEP) maps of chloroquine and chloroquine phosphate molecules.

Export to PPT

3.4 Molecular docking analysis

Molecular docking studies of chloroquine and chloroquine phosphate ligands were carried out with four structures of COVID-19 protein (PDB ID: 6 M03, 5R7Y, 5R81 and 6LU7). The two ligands were tested for drug-likeliness properties. Calculations were performed using the iGEMDOCK program through the generic evolutionary method (GA) and an empirical scoring function. Both ligands and target proteins structures were adapted with Discover Studio Visualizer software. All crystallographic water molecules were removed.

Our goal is to determine the modes of interaction of protein-ligand complexes while looking for favorable orientations for the binding of a ligand to a receptor (Duhovny et al., 2002; Seeliger and de Groot, 2010; Amin et al., 2010; Ahmed et al., 2013; Ghalla et al., 2018). In our case, the receptor represents the COVID-19 protein which has one or more specific active sites, more or less accessible. At each step, the interactions are affected and the best pose of the ligands is determined. 10 poses have been obtained; we have chosen the best pose with the lowest energy. These best poses, as presented in Fig. 7, were selected for investigating the different types of interactions that introduce a biological signal.

Close

Fig. 7

Orientation of chloroquine and chloroquine phosphate in the active sites of COVID-19 proteins.

Export to PPT

3.4.1

3.4.1 Chloroquine

The examination of Table 5 revealed that the chloroquine ligand presented the highest total energy score with the target protein 6 M03 which is equal to −81.866 kcal/mol. Note that the total energy is the sum of the three energies interactions: VDW, hydrogen band and electronic. Van der Waals interaction is a potential energy of attraction between two molecules. It represents the sum of the energies of Keesom, London and Debye. The H-bond represents an interaction between two electronegative atoms. Generally, the energy of an H-bond is of the order of a few tens of KJ/Mol. It varies between 1 and 60 KJ/mol for neutral fragments, and sometimes it can reach higher values for some covalent bonds. The last interaction is electronic; they always take very low values compared to the other two interactions.

Table 5 Docking results of chloroquine and chloroquine phosphate in COVID-19 protein.

Chloroquine
Ligands	6 M03	5R7Y	5R81	6LU7
Total energy	−81.866	−77.498	−68.514	−67.136
VDW	−75.581	−70.605	−65.014	−64.988
H-bond	−6.285	−6.893	−3.500	−2.147
Electronic	0	0	0	0
Affinity	−6.7	−6.6	−6.7	−6.1
Chloroquine phosphate
Ligands	5R7Y	6 M03	5R81	6LU7
Total energy	−99.119	−88.686	−84.817	−82.663
VDW	−66.409	−55.450	−79.862	−69.861
H-bond	−29.499	−30.505	−4.9547	−12.802
Electronic	−3.210	−2.731	0	0
Affinity	−4.5	−3.5	−3.5	−3.6

Chloroquine ligand posses the strongest van der Waals interaction E_VDW = -75.581 kcal/mol. The docking pose analysis showed that the chloroquine ligand is oriented with the VDW interactions surrounded by the chains of LEU-141, MET-165, PHE140, HIS163, GLN189, MET49, GLY143, THR25 and VAL42 binding residues in the 6 M03 protein. Also, it have the strongest H-bond interaction E_H-bond = -6.893 kcal/mol. The greater negative energy score suggests a more favorable binding mode. Table 6 presents the different interactions between the chloroquine ligand and proteins via the binding residues along with their bond length. Results obtained for protein targets show that the chloroquine ligand has bonded effectively with 6 M03 target sites with two remarkable carbon-hydrogen bond interactions. The mentioned compound is immensely bonded with active residues SER144 (Serine) and HIS164 (Histidine) by carbon-hydrogen bond interactions conduct to more antiviral activity. The first C—H bond interaction was identified between H₄₆ atom and SER144 binding residues and the distance was found to be 2.61 Å. The second C—H bond interaction was identified between H₂₇ and HIS164 with distance 2.27 Å. The hydrogen atom H₃₀ linked to HIS41 amino residues via an alkyl interaction with bond length equal to 4.11 Å. Also, Pi-Sulfur, Pi-Alkyl and Pi-Anion interactions were observed surrounded by the amino acids CYS145, LEU27 and GLU166, having distances 3.99, 4.28 and 4.55 Å, respectively. These results have been well described in Figs. 8 and 9. Furthermore, chloroquine molecule showed total energy score of −77.498 kcal/mol against 5R7Y protein with VDW interaction (−70.605 kcal/mol) and hydrogen bond energy (−6.893 kcal/mol). Regarding the two other proteins (5R81 and 6LU7), the interaction energies are slightly weaker in comparison with the other ligands but as even remain important. The docking calculations led to the following results: the total energies scores are equal to −68.514 kcal/mol and −67.136 kcal/mol for 5R81 and 6LU7, respectively. The van der Waals interactions were found to be E_VDW (for 5R81) = −65.014 kcal/mol and E_VDW (for 6LU7) = −64.988 kcal/mol. Additionally, the hydrogen bond interactions exhibiting values of −3.500 and −2.147 kcal/mol for 5R81 and 6LU7 receptors. In the chloroquine-5R7Y complex, a Pi-Anion and Pi-Sulfur interactions wrapped by the amino acids GLU166 and CYS145 were formed with bond lengths 4.42 and 4.03 Å. C₁₅ atom made two Alkyl interactions with A:CYS145 and A:LEU27 residues and having distances 3.99 and 4.07 Å. Also, C₁₅ interact with A:HIS41 via a Pi-Alkyl interaction (bond length = 3.88 Å). A:SER144 and A:HIS164 amino residues form two carbon-hydrogen bond interactions with H₄₆ and H₂₇ atoms. Their bonding distances are found to be 2.53 Å and 1.98 Å, respectively. In 5R81virus, A:MET165 and A:MET49 amino residues are involved in the alkyl interaction with C₁₀ and C₁₅ atoms having bond length 4.43 and 3.96 Å. Pyridine group formed Pi-Alkyl, Pi-Sulfur and Pi-Donor hydrogen bond interactions with A:LEU27 (5.13 Å), A:CYS145 (4.08 Å) and A:CYS143 (3.80 Å) residues, respectively. Another Pi-Alkyl interaction is also seen which contributed by A:HIS41 with C₁₅ atom, indicating distance 4.25 Å. For the last ligand 6LU7, the LEU141 (2.38 Å), the ASN142 (3.02 Å) and the HIS163 (2.47 Å) amino acids formed a C—H bond interactions with H₂₉, H₂₇ and H₂₈ atoms of chloroquine. In addition to these weak interactions there are two alkyl interactions; one between PRO168 residues and the Cl atom and the second one is in between CYS145 and the N₂ atom, indicating bond distance 3.63 and 4.35 Å, respectively. Subsequently, the H₃₀ atom exhibit a conventional-H bond interaction with GLU166 residues and bonding distance is 2.22 Å.

Table 6 Amino acid residues-chloroquine interactions.

Ligand	Target protein	Binding residue	Type	Atoms	Bond length (Å)	Interactions
Chloroquine	5R7Y	A:GLU166 A:CYS145 A:CYS145 A:LEU27 A:HIS41 A:SER144 A:HIS164	GlutamicAcid Cysteine Cysteine Leucine Histidine Serine Histidine	Benzene Pyridine C15 C15 C15 H46 H27	4.42 4.03 3.99 4.07 3.88 2.53 1.98	Pi-Anion Pi-Sulfur Alkyl Alkyl Pi-Alkyl Carbon-H bond Carbon-H bond
	6 M03	A:CYS145 A:GLU166 A:HIS41 A:LEU27 A:SER144 A:HIS164	Cysteine GlutamicAcid Histidine Histidine Serine Histidine	Pyridine Pyridine H30 H30 H46 H27	3.99 4.55 4.11 4.28 2.61 2.27	Pi-Sulfur Pi-Anion Alkyl Pi-Alkyl Carbon-hydrogen bond Carbon-hydrogen bond
	6LU7	A:LEU141 A:ASN142 A:HIS163 A:PRO168 A:CYS145 A:GLU166	Leucine Asparagine Histidine Proline Cysteine GlutamicAcid	H29 H27 H28 Cl N2 H30	2.38 3.02 2.47 3.63 4.35 2.22	C—H bond C—H bond C—H bond Alkyl Alkyl Conventional H-bond
	5R81	A:MET165 A:MET49 A:HIS41 A:LEU27 A:CYS145 A:CYS143	Methionine Methionine Histidine Histidine Cysteine GlutamicAcid	C10 C15 C15 Pyridine Pyridine Pyridine	4.43 3.96 4.25 5.13 4.08 3.80	Alkyl Alkyl Pi-Alkyl Pi-Alkyl Pi-Sulfur Pi-Donor H-bond

Close

Fig. 8

2D visual representations of chloroquine ligand-COVID-19 proteins.

Export to PPT

Close

Fig. 9

Different interactions between ligand and their receptor.

Export to PPT

Close

Fig. 9

Different interactions between ligand and their receptor.

Export to PPT

In order to upgrade the recognition of the interactions existing between receptor and ligand, the affinities of these complexes were calculated by using AutoDockTools (ADT) (Morris et al., 2008). These affinities describe the strength of a non-covalent interaction between the ligand and its target which binding to a site on its surface. It is premised on the numeral and the nature of the physicochemical interactions. As illustrated in Table 5; the affinities values (in ultimate value) of chloroquine are found to be in the order of 6.7 > 6.6 > 6.1 kcal/mol for (6 M03 and 5R81), 5R7Y and 6LU7, respectively.

3.5 Hybridization effect

Of course, each compound has its own characteristics that distinguish it from the rest. The chloroquine phosphate is initially made up of chloroquine. Evidently, the adding of other atoms in the geometry of the chloroquine has an influence on their stability. The chloroquine compound becomes more stable when adding the phosphate groups since the global minimum energy decreases. Moreover, the smallest dipole moment was obtained for the chloroquine whereas the highest one was obtained for the chloroquine phosphate. This increase shows that the chloroquine is harder before adding the phosphate groups and also it promotes the formation of hydrogen bonds. We also find that by adding phosphate group the gap energy decreases, which involves a high reactivity for the chloroquine phosphate. This decrease in gap energy makes the flow of electrons easier, so the molecule becomes soft and more reactive.