New In Silico Insights into the Application of (Hydroxy)Chloroquine with Macrolide Antibiotic Co-Crystals against the SARS-CoV-2 Virus

Abstract

In this in silico study, different pharmaceutical co-crystals based on (hydroxy)chloroquine with macrolide antibiotics (azithromycin, clarithromycin, or erythromycin A) were analyzed for the first time. These findings present a new molecular perspective and therefore suggest that the combination of (hydroxy)chloroquine/azithromycin, in the stoichiometric ratio of 1:1, as model co-crystal systems has less toxicity and is the most effective for inhibiting the new severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).

1. Introduction

Recently, enormous efforts have widely focused on facing a new coronavirus, the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a disease designated as COVID-19 [1]. The World Health Organization (WHO) declared the COVID-19 outbreak as a new pandemic on 11 March 2020 [2]. In order to contribute to the front line of the fight against this disease, researchers worldwide are moving to an investigation of an effective treatment regimen (based mainly on known drugs) against infection caused by the SARS-CoV-2 virus [1,2,3,4,5]. Most recently, the combination of drugs, for example, (hydroxy)chloroquine with azithromycin (AZ), has been used as a strategy and has also been clinically tested to treat the SARS-CoV-2 outbreak [6,7,8,9,10,11]. The possible benefits of this therapeutic option are still being debated [8,9,10,11]. It is additionally expected, however, that computational screening can aid in accelerating the discovery as well as the development of new and more effective drugs against SARS-CoV-2 [1]. To the best of our knowledge, no study reports the benefits of using pharmaceutical co-crystals to treat patients with infection caused by the SARS-CoV-2 virus. However, a simple consultation of the literature suggests that pharmaceutical co-crystal research has become a prevalent theme in the last few years [12,13,14,15,16,17,18]. In particular, pharmaceutic co-crystals are usually defined as multicomponent crystals that incorporate two or more drugs [13,14,18]. These co-crystals exhibit different and improved physical–chemical properties, and thus, they are widely considered good candidates for diverse therapeutic applications [14]. Further, a large variety of synthetic strategies has been successfully developed to obtain controlled pharmaceutical co-crystals with different structures. Notable examples include solid-state grinding, solution- and melt-crystallization, solvent evaporation, and so on [12,14,15,16,17]. Therefore, from this perspective, we believe that a computational model of co-crystal systems is essential to gain further insight into their biological activity, aiding directly in the emergence of new experimental research to confront COVID-19.

In this paper, we have selected (hydroxy)chloroquine (denoted as HQ and CQ, respectively) with macrolide antibiotics (such as AZ, clarithromycin (CL), or erythromycin A (ER)) in the stoichiometric ratio of 1:1, as a co-crystal model system for a computational screening study. First of all, these compounds were selected because of some preliminary clinical studies [6,7,8,9,10,11]. However, more recent results, which were released after the conclusion of this research, suggest the ineffectiveness of the use of CQ and HCQ in the treatment and prevention of COVID-19 [19]. Despite this fact, it is also important to emphasize that information on the mechanism of action of co-crystals as a new therapeutical option, especially from a computational point of view, is of course scarce in the literature. This is probably related to the difficulty of structurally treating these compounds in the active site of the target biomolecule. In addition, in this case, we observe that this study also represents advances in the computational treatment of such systems. Hence, this work can reveal the interaction profile and affinity of these co-crystals at short and long ranges, serving as a guide for future studies in this field. Second, AZ and CL are well-known semi-synthetic derivatives of ER, the first and most-known macrolide antibiotic to be isolated in 1952 [20,21,22]. These macrolide antibiotics are mainly used to treat patients with certain infections of the respiratory tract [20]. Further, it is well-known that the combination of AZ with antimalarial drugs (e.g., (hydroxy)chloroquine), shows some crucial benefits for the treatment of malaria [23]. A third relevant aspect of this study is that both (hydroxy)chloroquine and macrolide antibiotics usually crystallize in a monoclinic structure according to the CSD (Cambridge Structural Database [24], which is fundamental for the construction of these co-crystal model systems. Therefore, from this perspective, we have developed the first theoretical model of co-crystals applied in the treatment of SARS-CoV-2. We believe that this strategy can bring additional benefits to the use of this therapy as well as contributing to a better understanding of a molecular point of view.

Finally, our molecular target, the SARS-CoV-2 main proteinase (M^pro), is a key enzyme responsible for the processing of viral polyproteins [25,26,27]. Together with the papain-like proteases, M^pro is essential for the processing of polyproteins translated from viral RNA [28,29]. In a structural analysis, the M^pro enzyme consists of three domains. Located in a cleft between domains I and II, the M^pro active site presents a catalytic dyad formed by the conserved Cys145 and His41 amino acid residues. CQ and HCQ are well-known due to their antimalarial activity, but have possible effectiveness against COVID-19. HCQ inhibits terminal glycosylation of ACE2, the receptor for cell entry. ACE2 in the non-glycosylated state may less efficiently interact with the SARS-CoV-2 spike protein, making the viral entry difficult [30]. In this work, we evaluated the antimalarials CQ and HCQ in combination with other macrolide antibiotics, in the single crystal and co-crystal formats. This study can allow us to have a deeper understanding of the affinity of these co-crystals for the molecular target, in comparison with single crystals of antibiotics.

2. Computational Details

Six co-crystal model systems were prepared in this study, as shown in Figure 1, based on the combination of HCQ and CQ with the AZ, CL or ER antibiotics in the stoichiometric ratio of 1:1. As a strategy, the unit cell of these compounds was replicated for the investigated set of single and co-crystals, that is, their structural parameters were maintained as similar as possible. All of these model systems were then optimized (as the starting point for docking studies), and their structural and vibrational properties were fully evaluated at the PM6 theoretical level, using the Gaussian 09 package [31]. In order to better describe the electronic/optical properties of these model systems, we also employed the Kohn–Sham time-dependent density functional theory (TDDFT) using the B3LYP functional with the 6-31+G(d,p) basis set. In addition to the partial charges of the atoms, we calculated the natural bond orbitals (NBOs).

In the docking studies, the affinity of the single and co-crystals in the M^pro active site was investigated. The compounds were then docked inside the crystallographic structure of viral protein M^pro (PDB code 6LU7; resolution = 2.16 Å) [32], using the Molegro Virtual Docker program (MVD^®) [33], taking into account the same procedures employed previously [34,35,36]. According to our calculation protocol, the radius was about 20 Å, where the residues of the catalytic triad were kept flexible. Due to the nature of the docking methods, the calculations were carried out, generating approximately 50 poses (such as conformation and orientation) for each ligand studied.

In the MVD program, the MolDock score algorithm method used as a scoring function is based on the piecewise linear potential, which is fundamentally a simplified potential whose parameters are in turn fitted to protein–ligand structures, binding data scoring functions, and further extended in the Generic Evolutionary Method for molecular docking, including a new hydrogen bonding term as well as new charge schemes [33]. Along this line, the docking scoring function values, E_score, are usually defined by Equation (1):

$$E_{score}=E_{inter}+E_{intra}$$

(1)

wherein:

$$E_{inter}={\displaystyle \sum }_{i\epsilon ligand}{\displaystyle \sum }_{j\epsilon protein}\left[E_{PLP}\left(r_{ij}\right)+332.0 \frac{qiqj}{4r_{ij}^2} \right]$$

(2)

Note that E_PLP represents ‘‘piecewise linear potential’’, which consists of two different parameter sets, as described: one for the approximation of the steric term (i.e., Van der Waals) among atoms, as well as the other potential for hydrogen bonding. As can be seen, the second term is, of course, related to the electrostatic interactions among overloaded atoms. Typically, it is a Coulomb potential with a dielectric constant dependent on the distance (which can be approximately described as D(r) = 4r). Hence, the numerical value of 332.0 is then responsible for the electrostatic energy unit to be given in kilocalories per molecule [33].

E_intra is defined as the internal energy of each ligand. That is:

$$E_{intra}={\displaystyle \sum }_{i\epsilon ligand}{\displaystyle \sum }_{j\epsilon ligand}{E}_{PLP}\left(r_{ij}\right)+{\displaystyle \sum }_{flexiblebonds}A\left[1-\cos \left(m.\theta -{\theta }_0\right)\right] +E_{clash}$$

(3)

Note that the first part of the equation (double summation) is among all pairs of atoms in the ligand, except those connected by two bonds. Thus, in this equation, the second term denotes the torsional energy, where θ is the bond’s torsional angle. Hence, if several torsions could be determined, then each torsional energy is considered an average among them. The last term of this equation, E_clash, has (not taking into account infeasible ligand conformations) a penalty incurred of about 1.000 for those exhibiting a distance between two heavy atoms (e.g., with more than two bonds apart) smaller than 2.0 Å [33]. Thus, in the MVD program, the docking search algorithm is based on the whole the interactive optimization technique (inspired by Darwinian evolution theory), which implies a new hybrid search algorithm conveniently called guided differential evolution. As such, this hybrid combines a differential evolution optimization technique with a cavity prediction algorithm during the search process. As a result, this approach provides a simple, fast, and accurate description of potential binding modes (poses) [33,37,38].

Additionally, taking into account the experimental inexistence of acute toxicity data for the compounds employed in this study and in order to estimate the theoretical values of the lethal dose (denoted as LD₅₀), a rat model-based admetSAR predictor was used. This admetSAR approach is freely available online at http://lmmd.ecust.edu.cn:8000/ (accessed on 15 August 2020).

Based on our docking and ADMET results, we carried out molecular dynamics (MD) simulation using the GROMOS54A7 all-atom force field [39] using GROMACS 5.1 software [39,40]. The co-crystal-M^pro complex was then inserted into a 12 Å water box with the SPC solvation model, and sodium and chlorine ions were added as strategy for charge neutralization under periodic boundary conditions. In this study, the long-range electrostatic interactions were calculated assuming the Particle Mesh Ewald method with a cut-off of 12 Å and a time step of 1 fs. Initially, the system was minimized over 5000 cycles using the steepest descent algorithm and was then submitted to a production simulation of 60 ns at 300 K and 1 bar pressure. Finally, all trajectories were later analyzed using Visual Molecular Dynamics (VMD, version 1.9.3) [41].

3. Results and Discussion

We started this study by analyzing frequency calculations for these co-crystals and the thermochemical parameters obtained (at 298.15 K and 1.00 atm) from them using the standard relations in the gas phase, given in

Table 1. Analysis of thermodynamic parameters computed for the co-crystals at the PM6 theoretical level.

Co-Crystal	ΔE (kcal mol−1)	H° (kcal mol−1)	S° cal (mol K)−1	CV cal (mol K)−1
CQ/AZ	97.2	932.4	459.3	333.9
CQ/ER	48.9	929.2	462.3	336.1
CQ/CL	52.0	946.8	476.3	342.1
HQ/AZ	52.3	935.2	469.3	337.5
HQ/ER	0.0	931.6	471.5	340.2
HQ/CL	7.6	949.5	482.5	345.6

. As such, the comparison of the difference between the total electronic energy (ΔE), used to predict the stability of the co-crystals designed, presented in Table 1, suggests that the HQ/ER is the co-crystal that is more stable. In both HQ/AZ and CQ/AZ co-crystals, it is the least stable in their respective series. Thus, the stability order obtained for the computed co-crystals is HQ/ER > HQ/CL > CQ/ER > CQ/CL > HQ/AZ > CQ/AZ, i.e., suggesting that HQ-based co-crystals are more stable than their CQ counterparts (see Table 1).

The typical calculated infrared (IR) spectrum is shown in Figure 2a–f. According to these results, we observed only positive modes, suggesting that the structures proposed in this study are well optimized and therefore represent a minimum energy. Further, the IR-active modes confirm that co-crystals were obtained with a monoclinic-like structure. These results are consistent with experimental observations for both (hydroxy)chloroquine and macrolide antibiotics [1,12,42,43]. As can be seen in Figure 3, the computed UV-vis absorption spectrum of these co-crystals revealed only one well-defined peak. As such, the more intense UV band located in the range of 617 nm to 771 nm for the HQ-based co-crystals (AZ to ER) is due to a predominant HOMO-to-LUMO+2 transition. In comparison, for the CQ-based co-crystals (AZ to ER) this UV band was located from 624 nm to 955 nm (which is due to a predominant HOMO-to-LUMO+3 transition for CQ/AZ, and both CQ/CL and CQ/ER was best assigned to HOMO-to-LUMO+2 transition).

In all cases, we observed a red-shift UV-vis absorption spectrum of HQ-based and CQ-based co-crystals in the function of increasing structural polarization (summarized by the calculated dipole moment, as shown in

Table 2. Comparison of the dipole moment, HOMO–LUMO energy gap, exciton binding energy, hardness, softness, electronegativity, and electrophilicity values for the computed co-crystals at the TDDFT theoretical level.

Co-Crystal	Eg (eV)	EB (eV)	ECT	µ (D)	η (eV)	S (eV−1)	x (eV)	ω
CQ/AZ	1.22	0.26	0.009	−1.03	0.61	0.82	1.03	0.43
CQ/ER	1.26	0.07	0.141	−1.02	0.63	0.79	1.02	0.44
CQ/CL	1.46	0.47	0.337	−0.94	0.73	0.68	0.94	0.52
HQ/AZ	1.44	0.26	0.019	−1.10	0.72	0.69	1.10	0.46
HQ/ER	1.44	0.26	0.632	−1.13	0.72	0.70	1.14	0.45
HQ/CL	1.45	0.28	0.350	−1.12	0.73	0.69	1.12	0.46

). This was the main reason for the slight variation in the angles and bond lengths observed in the studied models (Figure 1).

Additionally, the calculated value of the HOMO–LUMO energy gap (E_g), exciton binding energy (E_B), hardness (η), softness (S), electronegativity (x), and electrophilicity (ω) values for the computed co-crystals at the TDDFT theoretical level are given in Table 2. According to the E_g values, we noticed an increase in the order CQ/AZ < CQ/ER < HQ/ER = HQ/AZ < HQ/CL < CQ/CL. Our results also suggest better chemical stability for the CQ/CL co-crystal. In addition, these HQ-based and CQ-based co-crystals have E_B values from 0.07 eV to 0.47 eV, which is typical of organic crystals containing aromatic groups [44].

Based on this analysis, supported by theory [45,46,47], we identified that HQ and CQ moieties are donor states, and the AZ, CL, or ER moieties are acceptor states. As such, the frontier molecular orbital analysis is illustrated in Figure 3 to view the energy and charge transfer process (denoted as E_CT) that occurs in these co-crystals. In general, this process has a profound impact on the functional properties of these model systems. The electronic excitation and natural bond orbital (NBO) analysis was performed using the TDDFT theoretical level and hence can be used to measure E_CT. Thus, in this perspective, the stabilization energy (E₂) associated with the delocalization of electrons between the electron donor NBO (i) and the electron acceptor NBO (j) was evaluated according to the following equation [48,49].

$$E_2=\Delta E\left(ij\right)=q_i \frac{{\left(F\left(i,j\right)\right)}^2}{{\epsilon }_j-{\epsilon }_i}=q_{i }\frac{F_{ij}^2}{\Delta \epsilon }$$

(4)

where q_i represents the orbital occupancy, ε_i and ε_j are diagonal elements (i.e., orbital energies), and F_ij is the off-diagonal NBO Fock matrix element. From these calculations, we calculated the E_CT for these co-crystals, as shown in Table 2. As expected, the HQ-based co-crystals had higher values of E_CT than their CQ counterparts. Consequently, in this case, we identified that the E_CT process was slower and hence can more easily induce symmetric structural polarization. For the CQ-based co-crystals, this result suggested a fast E_CT process that can more likely induce an asymmetric structural polarization.

According to the above co-crystal structure analysis, the molecular docking calculations were performed to adjust the ligands in the M^pro active cavity, evaluating the affinity among them. For this study, a cavity prediction algorithm based on a 3D box was employed in order to generate the binding sites. Thus, the volume of the calculated active cavity was approximately 130.56 ų.

In order to understand the interaction modes of these single and co-crystals within the SARS-CoV-2 M^pro active site, some parameters that contribute to affinity were investigated.

Table 3. Values of the parameters obtained by molecular docking calculations for the single and co-crystals within the SARS-CoV-2 M^pro active site through MVD^® software and toxicity prediction through ADMET.

Co-Crystal	Interaction Energy (kcal mol−1)	Toxicity (LD50)
AZ	−213.2	2.54 mol kg−1
CL	−205.0	2.72 mol kg−1
CQ/AZ	−275.3	2.74 mol kg−1
CQ/CL	−240.9	2.75 mol kg−1
CQ	−141.7	2.95 mol kg−1
CQ/ER	−229.8	2.75 mol kg−1
ER	−213.1	2.23 mol kg−1
HQ/AZ	−280.7	2.74 mol kg−1
HQ/CL	−259.6	2.75 mol kg−1
HQ	−162.6	2.63 mol kg−1
HQ/ER	−275.7	2.75 mol kg−1

shows the intermolecular interaction energy values and theoretical toxicity prediction for the ligands under study.

As shown in Table 3, all of the single and co-crystals studied interacted very well with the M^pro active site, with interaction energy values in the range of −141.7 to −280.7 kcal mol⁻¹. We also observed that the co-crystal formed by HQ and AZ had the highest stabilizing interaction energy (−280.7 kcal mol⁻¹), followed by the co-crystals HQ/ER (−275.7 kcal mol⁻¹) and CQ/AZ (−275.3 kcal mol⁻¹). The conformation of HQ/AZ in the M^pro active cavity is shown in Figure 4. Note that the single crystals of CQ and HQ had interaction energies of approximately −141.7 and −162.6 kcal mol⁻¹, respectively. However, when combined with AZ, these single crystals gained stability and had lower interaction energies. In fact, the single crystal of AZ was much more stable than those of CQ and HQ when docked in the active site, and their combination with AZ led to increased affinity for the viral enzyme. Regarding the co-crystals formed by the combination with AZ, we observed the following energy differences: HQ to HQ/AZ (118.1 kcal mol⁻¹) and CQ to CQ/AZ (133.6 kcal mol⁻¹). Based on our computations, CQ/AZ was the co-crystal that showed the more significant energy difference, followed by HQ/AZ, suggesting the suitability of combining AZ with CQ and HQ, thus boosting their therapeutic effects (Table S1).

According to Table 3, the combination of HQ with ER also led to stabilizing interaction energies, from −162.6 kcal mol⁻¹ for HQ to −275.7 kcal mol⁻¹ for HQ/ER. This trend was also observed for CQ; its interaction energy was stabilized to 88.1 kcal mol⁻¹ when combined with ER. Interestingly, the formation of the co-crystal HQ/ER was much more significant in energetic terms due to the larger energy difference between the single crystal and co-crystal, i.e., 113.1 kcal mol⁻¹. In addition, with respect to the co-crystals formed by the combination with CL, we observed the following energy differences: HQ to HQ/CL (97.0 kcal mol⁻¹) and CQ to CQ/CL (99.2 kcal mol⁻¹). In general terms, all macrolide antibiotics led to the formation of low-energy co-crystals, indicating that such co-crystals can be easily obtained experimentally.

With respect to the interactions observed in the molecular docking, we observed that most compounds (single and co-crystals) interacted with important residues in the site, such as the residues H41, C145, and E166. In the active site, H41 and C145 constituted the dyad catalytic of M^pro, being key residues for the enzyme’s inhibition process. Furthermore, E166 is also an important species for keeping the active conformation of the enzyme [3,50,51]. As a result, the interactions of all compounds investigated can be visualized in Figure S1 and Figure 5.

Usually, the efficacy of a drug can be attributed to many factors. Along this line, the absorption–distribution–metabolism–excretion–toxicity (ADMET) evaluations involve sequential and iterative assessments of the efficacy, pharmacokinetics, pharmacodynamics, metabolic, and toxicological properties in the model of therapeutic agents, contributing to the design of safer and more efficient drugs [52]. The ADMET analysis was carried out in order to evaluate the toxicity of our compounds, single and co-crystals, under investigation. According to our results in Table 3, the LD₅₀ value found for ER, 2.23 mol.kg⁻¹, indicates more toxicity in relation to other single and co-crystal compounds. However, we observed that its combination with CQ and HQ led to co-crystals of lower toxicity, as well as more stabilizing interaction energies, making these efficient species candidates for COVID-19 treatment. In addition, when using CL combined with CQ and HQ, there was a slightly decreased toxicity and significant interactions of lower energies in the M^pro active site. The most remarkable results were found for the co-crystals formed by the combination with AZ. The data reinforce our previous findings, which indicate that the co-crystals of CQ and HQ combined with AZ significantly improved the affinity and decreased toxicity, thus making them more potent therapeutic agents against SARS-CoV-2.

Additionally, in this study, the dynamic behavior of the HQ/ER-M^pro complex was investigated. The co-crystal HQ/ER had a stabilizing interaction energy and lower toxicity, according to the docking and ADMET results, respectively. This system was chosen for a deeper investigation through Molecular Dynamics (MD) simulation. Initially, we observed that the complex reached thermodynamic equilibrium after 45 ns of MD simulation (see Figure 6). The extracted frame, which was considered the representative conformational structure of the co-crystal throughout the MD simulation, corresponds to the average of the Root Mean Square deviation (RMSD) value. By analyzing the results of the RMSD graph, we noticed that most of the deviations from HQ/ER were larger during most of the MD simulation, suggesting that the co-crystal requires a period of time to reach a stable accommodation in the active site. By analyzing the graph of the hydrogen bond (Figure 6), we found that the co-crystal HQ/ER performed up to six hydrogen bond type interactions. However, more effective interactions were observed by the end of simulation (after 45 ns), after the co-crystal was well-accommodated and with stable interactions in the site. According to the total energy graph (Figure 6), the energy oscillated during all simulation periods; however, this oscillation was constant, without very large variations. In addition, based on the Root Mean Square Fluctuation (RMSF; Figure 6), we observed an oscillation of the backbone, protein, and sidechain, during the MD simulation. Based on the graph, this oscillation suggests conformational changes in the protein structure, allowing for a stable accommodation of the co-crystal in the site during the simulation.

It is important to emphasize that in silico research can contribute to a better physicochemical interpretation of a huge variety of biological processes at the molecular level, and in this direction, such strategies can contribute to accelerating the process of developing new drugs [53,54,55,56]. Usually, however, due to great structural complexity and high computational cost in the treatment of more realistic models, only a small part of the tar-get microorganism is used during these simulations. In other words, these results repre-sent only an initial effort of preclinical research. For this reason, only clinical trials can indeed confirm the effectiveness of such therapeutic options. Finally, from this perspective, we and other authors believe that in the future, the use of realistic models based on the complete structure of the target microorganism, as recently reported by Durrant et al. [57], will be dominant in this field of computational drug design. Undoubtedly, these advances will lead to a new era in computational drug design, providing better prediction power of how the structure will affect their biological activity under more realistic conditions, accelerating the discovery of new therapeutic options.

4. Conclusions

In summary, our theoretical findings are of great importance and provide a new molecular perspective on the effectiveness and toxicity of co-crystals formed between (hydroxy)chloroquine with macrolide antibiotics for the treatment of SARS-CoV-2 infection. Our findings also enable a new interpretation in for the intramolecular energy and charge transfer process in these co-crystals, which is critical to the rationalization of their functional properties, providing an innovative way to connect structural changes with electronic transfer kinetics. Overall, the co-crystals of CQ and HQ combined with AZ significantly improved the affinity and decreased toxicity, making them more potent therapeutic agents against SARS-CoV-2. However, it must be emphasized that the therapeutic potential benefits can only be determined in practice using rigorous clinical studies.