As a first step, we have performed DFT and TDDFT calculations for the nitro derivatives of compounds Q and QO to better understand of their electronic structure, spectroscopic properties, and chemical reactivity. Figure 3 shows the electrostatic surface potential for the optimised structure of all the nitro derivatives of Q and QO, which have screened in this in silico study. Hence, we can also see that the charge distribution mainly depends on the various orientations of the nitro groups - regions with negative potential (red) that act as an excellent electron acceptor - that were added to the Q and QO compounds, respectively, with the specific objective of conferring the most favorable interaction between the drug and the target. Note also that the nitro group increases the polarity of these compounds, which is an attractive characteristic for pharmacological applications.35 Additionally, nitro derivatives of QO compounds, in this case, have a more polarized structure. These slight structural changes are responsible for modulating the biological activity of these compounds, which may provide new clues for an in-depth interpretation of their microscopic behavior. These theoretical findings are consistent with the molecular docking simulations performed in this study.
All compounds were also identified in terms of computed IR-active modes and UV-vis absorbance spectroscopy, as we show in Figure S1. These results can easily be used to distinguish the isomers obtained. In parallel, a comparison of the difference between the total electronic energy (ΔE) for the computed Q and QO isomers, presented in Table S1, suggests that both N–4-Q and N–4-QO compounds in terms of energy are more stable. Additionally, the HOMO-LUMO gaps reveal a minor difference of 4.07 to 4.31 eV for nitro derivatives of Q and of 3.11 to 3.69 eV for nitro derivatives of QO, respectively (see Table S1). In this case, a lower HOMO-LUMO gap value for QO derivatives, in principle, suggests greater reactivity for these isomers compared to Q derivatives. Figure S2 shows the shape of molecular orbitals (MOs) for all ligands studied. A detailed analysis of composition and localization of the MO reveals that the HOMO energies are, in principle, insufficient to describe the chemical behavior of these ligands. From the frontier effective-for-reaction molecular orbital (FERMO) concept, the reactions that are driven by HOMO, and those that are not, can be better explained for such compounds.36–39 These findings are consistent with previous studies.38,39
To elucidate the modes through which our drug candidates interact SARS-CoV- 2, the crystal structure of the Mpro of the virus in complex with 6-(ethylamino)pyridine- 3-carbonitrile was downloaded for study.27 Once the enzyme had been prepared, the molecular docking protocol was initiated. In the first part of this investigation, re-docking calculations were performed using the MolAr software,29 with implementation of the AutoDock Vina program.28 To determine the ideal docking parameters, these re-docking calculations were performed according to the orientation and conformation adopted by the experimental co-crystallised active ligand present in the binding pocket. It is important to notice that the Mpro enzyme used in this work was found in its native form.
The small RMSD variation (0.94) obtained from the re-docking calculations, suggested that the program was able to correctly and efficiently simulate the experimental results for the respective ligands. This preliminary outcome indicated that the conformational deviation of the molecular docking technique was suitable for our purposes and that the method was highly sensitive and specific. The re-docking overlap is presented in Figure 1. To simulate the modes through which our drug candidates interact with the SARS-CoV–2 Mpro enzyme, we employed the best parameters provided by the data from the re-docking study carried out with the co-crystallised active ligand. All the computed interaction energy results are displayed in Table S2 and S3 in supplementary material.
As shown in Table S2, all the drug candidates studied (i.e., nitro derivative of Q and QO) interacted well with Mpro active site, with interaction energy values in the range of –4.3 to –5.0 kcal.mol–1, respectively. Some of the nitro-QO compounds, such as the inhibitors N–4-QO, N–9-QO, N–8-QO, together with QO, had slightly more stabilising interaction energy values than those of their corresponding nitro-quinolines (Table S2). In general, it is noteworthy that the compounds studied had a greater affinity for Mpro than the co-crystallised ligand did (the latter showing an interaction energy value of –3.9 kcal mol–1). In order to assess the potential of such findings, using the same protocol, docking procedures were performed with the commercial drugs CQ and HCQ, which are currently adopted for the treatment of SARS-CoV–2 infection, and their interaction energy values were found to be –2.8 and –2.3 kcal mol–1, respectively. A remarkable trend could be observed from these outcomes. Note that all of our drug candidates presented greater interaction energy values than CQ and HCQ, with a significant energy difference, of up to 2 kcal mol–1. Additionally, our study showed that the many of the nitro-QO compounds led to a more stabilising interaction energy in the Mpro active site. Based on these findings, we also investigated the chloroquine and hydroxychloroquine N-oxides forms (denoted as CQO and HCQO), which displayed a significant improvement in interaction energy values of –3.0 and –3.1 kcal mol–1, respectively. Interestingly, the interaction energy of HCQO was almost 1 kcal mol–1 more stabilizing than the HCQ. This trend was deeply analyzed using molecular dynamics (MD) techniques. The influence of the N-oxide group was also investigated at different sites and through different combinations for the CQO and HCQO compounds (Table S3). According to that table, with all combinations investigated, we can observe that no improvement in interaction energy was detected for CQO. On the other hand, for HCQO, the presence of the N-oxide group at some sites led to slightly more favorable interaction energies. See table S3 for more details. Herein, our main goal was to determine whether the inhibitors studied could target the Mpro enzyme. The molecular docking pose of each drug candidate indicated that they could indeed fit accurately within the substrate-binding pocket.
In the case of SARS-CoV–2 virus Mpro enzyme, the protomer is composed of three domains, as commented previously (see Figure 1). The enzyme has a Cys145–His45 catalytic dyad, and the substrate-binding pocket is known to be located in a cleft between domains I and II.40 Hence, the structural features determined from these data are important for guiding our assessment of the interaction modes of the inhibitors in the Mpro active site. As shown in Figure 4, the N–4-QO performed hydrogen bonds with all the residues and the water molecule of the catalytic triad. In fact, these specific interactions constitute one of the parameters analysed in this docking study. This same trend is not observed for inhibitor N–4-Q, suggesting that the N-oxide version of this ligand adopts a more favourable conformation which allows for its interaction with the catalytic triad, resulting in a slightly more stabilizing interaction energy. Similarly, the interactions performed by the other ligands can also be observed in Figure 4.
From the molecular docking calculations, it was possible to deduce that our drug candidates had a more stabilizing interaction energy effect than CQ and HCQ in the Mpro binding pocket. To better assess the interaction modes of our inhibitors, N–4-Q and N–4- QO were chosen as representatives of the set for MD simulations. Likewise, the same calculations were performed for CQ and HCQ and their N-oxides CQO and HCQO (see Figure S3).
Additionally, in this study, the dynamic behavior of complexes Mpro/N–4-Q, Mpro/N–4-QO, Mpro/CQ, Mpro/CQO, Mpro/HCQ, Mpro/HCQO inside the SARS-CoV–2 Mpro enzyme was investigated. Initially, we have observed that all complexes reached thermodynamic equilibrium during the 10 ns of MD simulation (see Figures S4-S9 in the supplementary material). The extracted frame, which was considered the representative conformational structure for all inhibitors throughout the MD simulation, corresponds to the average of the RMSD value. By analyzing the results of the RMSD graphs, it was observed that most of the deviations from the N–4-Q and N–4-QO structures were very small, not exceeding 0.5 Å, i.e., these ligands are well-accommodated in the SARS-CoV- 2 Mpro active site according to Figures S4 and S5.
To get more insights into the intrinsic reactivity of each these ligands, in this study, we have performed the analysis of the strain effect along the MD simulation. Since these factors have a pivotal role and affect the reactivity of these ligands.41,42 In the present study, the strain effect along the MD simulation can be clearly visualized by the overlap of the initial (blue) and representative (red) structures obtained after 10 ns of simulation, as shown in Figure 5. Based on Figure 5, we note that our compounds N–4-Q and N–4- QO showed a slight bending at the quinoline ring (strain), which makes this compound in principle more reactive, resulting in a small oscillation in according to the RMSD graphs (Figures S4 and S5). Importantly, this trend is essential because it indicates a low variation of strain (deformation of the ligand along simulation), reaching a stabilizing conformation more quickly. On the other hand, due to larger molecular mass and bulk of the CQ, HCQ and their corresponding oxides, there was a very higher variation of strain (Figures 5 (c—f)). We believe that this slight strain can induce a high intrinsic reactivity for these ligands.
As shown in Figure 6 (a), the N–4-Q compound performed hydrogen bond interactions with His41 (2.78 Å) and hydrophobic interactions with Leu27, Asn28, Pro39, Val42, Ile43, Cys44, Met49, Asn119, Gly143 and Gys145, respectively. These interactions are essential in inhibiting the enzymatic activity of Mpro and are in accordance with some other studies,43–45 as well as the docking results of this work. By analyzing the graph of hydrogen interactions (see Figure S4), we found that the compound N–4-Q performed up to four hydrogen bond type interactions. However, there was only one effective interaction that occurred during the entire 10 ns of MD simulation, which is according to the pharmacophoric map (see Figure 6 (a)). On the other hand, the N–4-QO compound is stabilized by two hydrogen bonds with His41 (3.17 Å) and His164 (2.72 Å) and also hydrophobic interactions with Leu27, His41, Cys44, Thr45, Ser46, Met49, Cys145, His164, Met165 and Glu166, respectively, as shown in Figure 6 (b). According to Zhang and coworkers,45 in the catalytic site, the residues Glu166, His41, and Gys145, respectively, are key species of the target protease. Thus, the interaction of these amino acids with inhibitors is essential for blocking the enzymatic activity of Mpro. Additionally, it is observed that the N–4-QO can make up to three bonds during the trajectory; however, occurs only one hydrogen bond in most of the entire simulation (Figure S5).
In the case of the dynamic behavior of both CQ and CQO compounds, we have observed that the structure remained stable during the 10 ns of simulation, as shown in Figures S6 and S7. In this case, the RMSD values found were around 2.5 Å. Although the RMSD deviation was high when compared to N–4-Q and N–4-QO compounds, they are coherent since the chemical structures of CQ and CQO are bulkier and had a higher molecular mass, as well as several rotational atoms. Consequently, there is a change in the conformation (Figure 6 (c) and (d)), further increasing the degrees of freedom of the inhibitors, and therefore is expected a more significant oscillation in the RMSD (see Figures S6 and S7). Through the pharmacophoric map, as shown in Figure 6 (c), it was observed that CQ performed hydrogen bond interactions with Ser46 (2.84 Å), Gln189 (2.09 Å) and Thr190 (2.67 Å) and had hydrophobic interactions with Thr25, Leu27, His41, Met49, Thr90, Ser144, Cys145, His164, Met165, Glu166, Leu167, Pro168, Ala191 and Gln192, respectively. It is essential to highlight that most of these residues are key species in catalytic activity.45 Note that the CQ can perform up to three Hydrogen bonds, with only one present in the entire simulation (see also Figure S6). On the other hand, in the case of CQO, this inhibitor performed several hydrophobic interactions, specifically with the residues His41, Val42, Ser46, Leu50, Cys145, His164, Met165, Glu166, Gln189 and Thr190 (see Figure 6 (d)). These results are in accordance with the Hydrogen bond graph, since, no hydrogen bonds are observed during the 10 ns of simulation (Figure S7).
Likewise, the HCQ compound remained stable at the Mpro active site only after 8.5 ns of simulation, mainly due to many conformational changes (Figure S8), such as in relation to the amine group that underwent rotation, resulting in a more energetically favourable conformation compared to its initial chemical structure, i.e., thus decreasing the RMSD value (Figure S8 and Figure 5 (e)). However, the final configuration of CQ had less strain than HCQ. Therefore, we can speculate that the minimum strain, which is associated with a low RMSD value, is directly related to the toxicity of these compounds.
As shown in Figure 6 (e), the HCQ exhibit one hydrogen bond with Tyr54 (2.89 Å) and several hydrophobic interactions with the residues His41, Cys44, Met49, Met165, Glu166, Leu167, Pro168, Phe181, Val186, Asp187, Gln189, Thr190 and Gln192, respectively. Through the Hydrogen bond graph, it was observed that for this compound after equilibration in the active site, up to two hydrogen bonds could be accomplished (Figure S8). While the HCQO compound showed the RMSD value of around 2.5Å, due to the structural distortions in the N-diethyl-pentane portion during the 10 ns of simulation (Figure S9). Hence, the most significant change happened approximately 6.3 ns, wherein the ligand underwent a 180° rotation, resulting in a less energetically favorable conformation compared to its initial chemical structure (Figure 5 (f)). Also, this compound performed three hydrogen bond interactions with His164 (2.61 Å) and Glu166 (1.82 Å) and hydrophobic interactions with Pro39, His41, Cys44, Ser46 Glu47, Met49, Cys145, Met165 and Gln189, respectively, as we shown in Figure 6 (f). Considering the hydrogen graph (Figure S9), the HCQO can make up to five hydrogen bonds.
From the MD simulation, we have estimated the interaction energies for all cases studied: N–4-Q (–13.65 kJ/mol), N–4-QO (–60.27 kJ/mol), CQ (–28.11 kJ/mol), CQO (-16.50 kJ/mol), HCQ (–54.88 kJ/mol), HCQO (–72.21 kJ/mol). An important outcome observed in this study is that the majority of the N-oxide compounds had an energetically more favorable affinity at the Mpro active site than their Q counterparts. These findings are consistent with the molecular docking calculations. The existence of intermolecular interactions strongly guides these trends. Note that the N–4-QO performed more hydrogen bonds than N–4-Q, as shown in Figure 6. This fact helps explain the more stabilizing interaction energy found for N–4-QO. From the pharmacophoric maps shown in Figure 6, the accomplishment of hydrogen bonds, along with the hydrophobic interactions, are key to understand the biological activity of these inhibitors.
Yet, the success of a novel drug candidate is commonly attributed to diverse factors, including their bioactivity, rich pharmacokinetic (PK) and pharmacodynamics (PD) profiles, as well as toxicity. It would be therefore of huge interest to investigate these properties in the preliminary stages to in silico design of safer and more efficient drugs. Hence, the ADMET evaluations involve sequential and iterative assessments of the efficacy, PK, PD, metabolic and toxicological properties in the model of potential drug candidates.46 From the ADMET results, the theoretical parameters of toxicity (LD50) were obtained for each compound: N–4-Q (2.53 mol/kg), N–4-QO (2.56 mol/kg), CQ (2.95 mol/kg), CQO (2.68 mol/kg), HCQ (2.66 mol/kg), HCQO (2.69 mol/kg). We can observe that the parameter toxicity slightly varied from N–4-Q to N–4-QO, suggesting that the toxicity of these compounds is essentially equal. Similarly, this trend also is observed for HCQ and its corresponding N-oxide (HCQO). On the other hand, we have noticed a more significant variation for CQ and CQO compounds, indicating that CQO theoretically presents a higher level of toxicity. In addition, these molecular calculations also showed that HCQ is more toxic than CQ. Yet, this trend does not corroborate with previous experimental results.47 It is essential to highlight that the molecular results obtained do not take into account the effects of the counterion and, for this reason, suggest a different trend to the experimental findings previously reported.47 In particular, this divergence most-likely is related to the fact that the commercially used HCQ is a salt- based on hydroxychloroquine sulfate, while the CQ used is a salt-based on chloroquine diphosphate. It additionally is well-known that the counterion has a substantial effect not only on its biological activity but also on the toxicity of such compounds as well.48–51 Therefore, our molecular results indicate that the presence of phosphate groups contributes to increasing the toxicity of CQ in the treatment of the SARS-CoV–2 infection. Consequently, we can conclude that for the same type of salt used, it is expected CQ to be less toxic than HCQ, according to the molecular trend observed in this study.