New in silico insights into the application of the (hydroxy)chloroquine with macrolide antibiotics co-crystals against the SARS-CoV-2 virus

In this in silico study, the different pharmaceutical co-crystals based on the (hydroxy)chloroquine with the macrolide antibiotics (azithromycin, clarithromycin, or erythromycin A) was analyzed for the rst time. These ndings 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-crystals systems have less toxicity as well as is the most effective for inhibiting the new coronavirus SARS-CoV-2.


Introduction
Recently, enormous efforts have widely been focused on facing a new coronavirus, SARS-CoV-2, disease designated as COVID-19. [1] The World Health Organization (WHO), declared COVID-19 outbreak as a new pandemic, on March 11, 2020. [2] In order to contribute to the front line of the ght against this disease, researchers worldwide are moving the 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 as, for example, the (hydroxy)chloroquine with azithromycin (AZ) has been used as a strategy and also clinically tested to treat the SARS-CoV-2 outbreak. [6][7][8][9][10][11] As yet, the possible bene ts 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 bene ts 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-crystals research has become a prevalent theme in the last years. [12][13][14][15][16][17][18] Particularly, pharmaceutic co-crystals are usually de ned as multicomponent crystals that incorporates 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] Also, 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, in this perspective, we believe that a computational model of co-crystals 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 the (hydroxy)chloroquine with the macrolide antibiotics (such as AZ, clarithromycin (CL), or erythromycin A (ER)), in the stoichiometric ratio of 1:1, as a co-crystal model system to a computational screening study. First of all, these compounds were selected because of some preliminary clinical studies. [6][7][8][9][10][11] Second, the AZ and CL are well-known semi-synthetic derivatives of ER, the rst and most-known macrolide antibiotic to be isolated in 1952. [19][20][21] These macrolide antibiotics are mainly used to treat patients with certain infections of the respiratory tract. [20] Also, it is well-known that the combination of AZ with antimalarial drugs (e.g., (hydroxy)chloroquine), shows some crucial bene ts for the treatment of malaria. [22] 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), [23] which is fundamental for the construction of these cocrystals model systems. Therefore, in this perspective, we have developed the rst theoretical model of co-crystals applied in the treatment of SARS-CoV-2. We believe that this strategy can bring additional bene ts to the use of this therapy as well as contributing to a better understanding of a molecular point of view.

Computational Details
Six co-crystal model systems were prepared in this study, as shown in Figure 1, and are based on the combination of (hydroxy)chloroquine (denoted as HQ and CQ, respectively) with the AZ, CL or ER antibiotics, in the stoichiometric ratio of 1:1. As a strategy, the unit cell of these compounds was here replicated to the investigated set of single-and co-crystals, that is, their structural parameters have maintained close as possible. All 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, from using the Gaussian 09 package. [24] In order to better describe the electronic/optical properties for these model systems have also been 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 was also elucidated by the natural bond orbital (NBO) calculations.
In the docking studies, the a nity 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 Å), [25] using the Molegro Virtual Docker program (MVD®), [26] taking into account the same procedures employed previously. [27][28][29] According to our calculation protocol, it was considered a radius of about 20 Å, where the residues of the catalytic triad were kept as exible. Due to the nature of the docking methods, the calculations were carried out, generating approximately 50 poses (hence 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 fundamentally is a simpli ed potential whose parameters are in turn tted to protein-ligand structures, binding data scoring functions and further extended in Generic Evolutionary Method for molecular docking, including a new hydrogen bonding term as well as new charge schemes. [26] Along this line, the docking scoring function values, E score , are usually de ned by Eq.

1:
Note that the E PLP represents ''piecewise linear potential'', which consists of two different parameter sets, as described forward: one for the approximation of the steric term (i.e., Van der Waals) among atoms, as well as the other potential for the 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, for this, the numerical value of 332.0 is then responsible for the electrostatic energy unit to be given in kilocalories per molecule, as well. [26] E intra is de ned as the internal energy of each ligand. That is: Note that the rst part of the equation (double summation) is among all pairs of atoms in the ligand, taking off 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 as an average among them. Being that the last term of this equation, E clash , attributes (not taking into account infeasible ligand conformations) a penalty incurred of about 1.000 in those exhibiting the distance between two heavy atoms (e.g., with more than two bonds apart) is smaller than 2.0 Å. [26] Thus, in the MVD program, the docking search algorithm is based on the whole the interactive optimization techniques (inspired by Darwinian evolution theory), which implies a new hybrid search algorithm conveniently so-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 leading that way a simple, fast, and accurate description of potential binding modes (poses). [26,30,31] 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 50 ), a rat model-based admetSAR predictor was used. Hence, this admetSAR approach is freely available online at http://lmmd.ecust.edu.cn:8000/.

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, and are given in Table 1. As such, the comparison of the difference between the total electronic energy (ΔE), used to predict the stability of these co-crystals designed, presented in Table 1, suggests that the HQ/ER is the co-crystal more stable. In both HQ/AZ and CQ/AZ co-crystals, we can evidence that 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 well more stable than their CQ counterparts (see Table 1). On the whole, the typical calculated Infrared (IR) spectrum is shown in Figure 2 (a-f). According to these results, we have observed only positive modes suggesting that the structures proposed in this study are well optimized and therefore represent a minimum of energy. Also, the IR-active modes con rm of obtaining co-crystals with a monoclinic-like structure. These results are consistent with experimental observations for both (hydroxy)chloroquine and macrolide antibiotics, respectively. [1,12,32,33] As can be seen in Figure 3, the computed UV-vis absorption spectrum of these co-crystals reveals only one wellde ned peak. As such, the more intense UV band located in the range of 617 nm to 771 nm for the HQbased 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 is located in the around of 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 is best assigned to HOMO-to-LUMO+2 transition), respectively.
On the other hand, in all cases, we have 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 calculated dipole moment, as shown in Table 2). 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 cocrystals at the TDDFT theoretical level are given in Table 2. According to the E g values, we notice 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 an E B values range of 0.07 eV to 0.47 eV, which is typical of the organic crystals containing aromatic groups. [34] Based on this analysis, supported by theory, [35][36][37] we identify 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 its functional properties of these model systems. Particularly the electronic excitation and natural bond orbital (NBO) analysis was performed using TDDFT theoretical level and hence can be used to measure E CT . Thus, in this perspective, the stabilisation energy (E 2 ) associated with the delocalisation of electrons between the electron donor NBO (i) and the electron acceptor NBO (j) is evaluated according to the following equation. [38,39] 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 have calculated the E CT for these co-crystals, as shown in Table 2. As expected, the HQ-based co-crystals have higher values of E CT than their CQ counterparts. Consequently, in this case, we identify that the E CT process is slower and hence can more easily induce symmetric structural polarization. For the CQ-based co-crystals, this result suggesting 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 a nity 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 Å 3 .
In order to understand the interaction modes of these single-and co-crystals within the SARS-CoV-2 M pro active site, in particular, some parameters that contribute to a nity were investigated. Table 3 shows the intermolecular interaction energy values and theoretical toxicity prediction for the ligands under study. As shown in Table 3, all the single-and co-crystals studied interacted very well with M pro active site, with interaction energy values in the range of -141.7 to -280.7 kcal mol -1 . We can also observe that the cocrystal formed by unities of HQ and AZ showed the most stabilizing interaction energy (-280.7 kcal mol -1 ), followed by the co-crystals HQ/ER (-275.7 kcal mol -1 ) and also CQ/AZ (-275.3 kcal mol -1 ). 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 revealed interaction energies of approximately -141.7 and -162.6 kcal mol -1 , respectively. Yet, when combined with AZ, these single-crystals gain stability and lower interaction energies. In fact, the single-crystal of AZ is much more stabilized than those of CQ and HQ when docked in the active site, and their combination with AZ leads to increased a nity toward the viral enzyme. Regarding the co-crystals formed by the combination with AZ, we have the following energy differences: HQ to HQ/AZ (118.1 kcal mol -1 ) and CQ to CQ/AZ (133.6 kcal mol -1 ). Based on our computations, CQ/AZ was the co-crystal that showed the more signi cant energy difference, followed by HQ/AZ, suggesting the suitability of combing AZ with CQ and HQ, thus boosting their therapeutic effects.
According to Table 3, note that the combination of HQ with ER also led to quite stabilizing interaction energies, from -162.6 kcal mol -1 for HQ to -275.7 kcal mol -1 for HQ/ER. This trend is also observed for CQ, in which its interaction energy was stabilized by 88.1 kcal mol -1 when combined with ER. Interestingly, the formation of the co-crystal HQ/ER was much more signi cant in energetic terms, due to the larger energy difference between single-crystal and co-crystal, of 113.1 kcal mol -1 . In addition, with respect to the cocrystals formed by the combination with CL, we have the following energy differences: HQ to HQ/CL (97.0 kcal mol -1 ) and CQ to CQ/CL (99.2 kcal mol -1 ). In general terms, all macrolide antibiotics led to the formation of low energy co-crystals and, therefore, support that such co-crystals can be easily obtained experimentally.
With respect to the interactions observed in the molecular docking, we observe 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 constitute 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,40,41] As a result, the interactions of all compounds investigated can be visualized in Figure S1 and Figure 5.
Usually, the e cacy of a drug can be attributed to many factors. Along this line, the absorptiondistribution-metabolism-excretion-toxicity (ADMET) evaluations involve sequential and iterative assessments of the e cacy, pharmacokinetics, pharmacodynamics, metabolic and toxicological properties in the model of therapeutic agents, and hence contributing to design safer and more e cient drugs. [42] The ADMET analysis was carried out in order to evaluate the toxicity of our compounds, singleand co-crystals, under investigation. According to our results in Table 3, the LD 50 value found for ER, 2.23 mol.kg -1 , indicates more toxicity in relation to other single-and co-crystal compounds. However, we can observe that its combination with CQ and HQ leads to co-crystals of lower toxicity, as well as more stabilizing interaction energies, making these e cient species candidates for the COVID-19 treatment. In addition, when using CL combined with CQ and HQ, there was a slightly decreased toxicity and signi cant interactions of lower energies in the M pro active site. The most remarkable results were found for the cocrystals formed by the combination with AZ. The data reinforce our previous ndings, which indicate that the co-crystals of CQ and HQ combined with AZ signi cantly improved the a nity and decreased toxicity, thus making them more potent therapeutic agents against SARS-CoV-2.

Conclusions
In summary, our theoretical ndings are of large importance and provided a new molecular perspective about the effectiveness and toxicity of co-crystals formed between (hydroxy)chloroquine with the macrolide antibiotics for the treatment of SARS-CoV-2 infection. Our ndings also enable a new interpretation in-depth for the intramolecular energy and charge transfer process in these co-crystals, which is critical to the rationalization of their functional properties, it provides an innovative way to connects structural changes with electronic transfer kinetics. Overall, the co-crystals of CQ and HQ combined with AZ signi cantly improved the a nity and decreased toxicity, making them more potent therapeutic agents against SARS-CoV-2. Yet, it must be emphasized that this therapeutic option potential bene ts can only be determined in the practice from rigorous clinical studies.