Comparative Computational Analysis of Dirithromycin and Azithromycin in Search for a Potent Drug against COVID-19 caused by SARS-CoV-2: Evidence from molecular docking and dynamic simulation

Department of Chemistry, College of Science and Arts in Baljurashi, Albaha University, Albaha, Saudi Arabia Faculty of Science of Sfax, Department of Chemistry, University of Sfax, B.P. 1171, Sfax 3000, Tunisia University of Manouba, ISBST, BVBGR-LR11ES31, Biotechpole Sidi Thabet, 2020 Ariana, Tunisia Department of Biology, University of Hail, College of Science, P.O. Box 2440, Ha’il 2440, Saudi Arabia Laboratory of Bioressources: Integrative Biology and Recovery, High Institute of Biotechnology-University of Monastir, Monastir 5000, Tunisia Department of Chemistry, College of Science, Qassim University, Buraidah 51452, Saudi Arabia University of Monastir, Faculty of Science of Monastir, Department of Chemistry, Monastir 5019, Tunisia Department of Pharmacy, University of Salerno, Via Giovanni Paolo II, 132, Fisciano, Salerno 84084, Italy Laboratory of Genetics, Biodiversity and Valorisation of Bioressources, High Institute of Biotechnology-University of Monastir, Monastir 5000, Tunisia


Introduction
The recent emergence of fatal coronavirus disease 2019  in Wuhan, China caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been actively spreading worldwide and becomes a serious threat to human health infecting about 193 million people as of 22 July 2021 (1,2). It has been reported that SARS-CoV-2 enter cells via binding to angiotensin-converting enzyme 2 (ACE2), mediated via the viral surface spike glycoprotein (S protein) followed by its priming by the transmembrane protease serine 2 TMPRSS2 (3). Together with blocking both TMPRSS2 and the endosomal cysteine proteases cathepsins, B and L activity inhibits entry of SARS-CoV-2 in vitro because they constitute the main targets for antiviral intervention (4). These therapeutic targets may act only on the initial phase of SARS-CoV-2 due to the release of chronic inflammation phase in the latter phase. Considering the high mortality rate of COVID-19 (2.15% all over the world), and despite the impressive efforts that have been applied, no approved prophylactic or therapeutic agent except some vaccines which not been approved to combat the severity of the new mutant strains of SARS-CoV-2. Besides vaccine development, the current pandemic demands an urgent call for the discovery of new potential molecules. In this context, despite the use of modern medicine to natural product and their bioactive molecules which possess several pharmacological properties with favorable efficacy and tolerable toxicity (5)(6)(7)(8)(9)(10)(11). Also, macrolides remain a major option to develop new drugs. This has prompted the world to respond with the development and discovery of potential inhibitors like macrolides with highly desirable anti-COVID-19 (12,13). Macrolides are a group of bacteriostatic antibiotics composed of 12-to 16-atom large lactone rings containing one or more deoxy sugars with usually attached cladinose and desosamine. Lastly, the antiviral effects of macrolides have attracted much attention because of their broad spectrum of activity displaying anti-inflammatory and antibacterial potency, commonly associated with respiratory tract infections in patients and viral pneumonia related to influenza and other viruses (14). This is due probably to their immune-modulatory as well as antiinflammatory actions by variably affecting cytokine expression, promoting the induction of immunoglobulin antibodies, and therefore reducing the susceptibility of respiratory viral infections (15). Among them, rhinoviruses, respiratory syncytial virus, and influenza virus have been reported to be inhibited by clarithromycin and azithromycin (14,16). Also, azithromycin can inhibit Zika and Ebola viruses as well as can reduced viral load and can induce epithelial interferon expression that contributes to the prophylactic effect of this drug in reducing exacerbations of chronic obstructive pulmonary (17)(18)(19). Influenza progeny virus replication was noticeably inhibited by azithromycin before infection (16). On the other hand, azithromycin, clarithromycin, erythromycin, fidaxomicin, and telithromycin were proved to be effective against influenza complications of respiratory viral infections (20)(21)(22)(23). Recent studies have shown that patients treated with azithromycin in combination with hydroxychloroquine have been demonstrated to exhibited a virological cure for . Thus, due to their above properties, and since the safety profile of the selected macrolide has already been approved, their repurposing is the best approach to find therapeutics anti-SARS-CoV-2.
To the authors' knowledge, no study has been made on the targeting of the macrolide, dirithromycin for COVID-19. Also, this is the first report that has taken account of SARS-CoV-2-15 receptors. For this, molecular docking and dynamic simulation were performed to screen the macrolide, dirithromycin against SARS-CoV-2-15 receptors to elucidate the potential inhibitors of its catalytic domain and the results were compared to the well-known azithromycin. These findings may be beneficial for the controlling of the propagation of the Novel Coronavirus  as well as the prediction of a new medication.

Materials and Methods Ligands design and docking
Dirithromycin and azithromycin were studied in the present work ( Figure 1). The 3D Structure of TMPRSS-2 required for the docking in this work was done using the I-tasser software (25) the best model with the highest C-score was used. Different methods were used for the design of the 3D structures of the targeted molecules required for docking in this work. Five small peptides were assessed by the PEPstrMOD software (26)(27)(28) and refined using Chimera software 1.14 using the force filed AMBER ff14SB. The other eighteen compounds were designed following the conversion of the 2 D coordinates of each molecule available on PubChem (29), made with the highest degree of precision available with the PyMOL tool (30). Docking between fifteen receptor and twenty-three molecules was carried out by using the AutoDock 1.5.7rc1 software (31,32). The coordinate of each receptor, free or complexed with inhibitor was extracted from his crystal structure available in the protein data bank (33) During the docking procedure, small peptides and organic compounds were used as ligands, only amide bonds were defined rotatable and almost all other bonds were defined as no rotatable. All receptors were kept rigid. Grid maps representing target proteins were constructed with different dimensions depending on the active site of the target protein ( Table 1).
The fourteen free three-dimensional (3D) molecular structures of the receptors or complexed with inhibitor and the TMPRSS-2 were visualized using the molecular visualization software PyMOL (30,34) and the ligand-protein complex interaction obtained was visualized by the Python Molecular viewer 1.5.6 (24) and the two-dimensional representation of the twodimensional representation was performed by Discovery Studio Visualizer 20.1.0 (35,36). Molecular dynamic simulations based on the docking results were applied to the azithromycin-furin and dirithromycin-furin complexes, and to the protein alone using NAMD v.2.13 in the CHARMM36 force field. The dynamic study was used to (i) determine the consistency and stability of the obtained complexes, (ii) follow the conformational changes of the pre-set simulation condition and (iii) judge their efficiency to bind and inhibit the target protein. The charmmgui.org module was used to construct the topology and the parameter files for the two inhibitors such us azithromycin and dirithromycin, used as input into the VMD software (37) to generate the .pdb coordinate files and the .psf structure files. The "add solvation box" module in the VMD modeling menu was used to generate the structure and coordinate files of the solvation structure required to run the molecular dynamics simulation. The solvation model generated by the VMD software used is a water cube covering the total protein-free or the whole complex and our dynamics simulation method is composed of three steps. The first is a minimization step performed for 1000 cycles to stabilize the complexes. The last frame of this step was used to start the second step of DM simulations for 140 ps that correspond to 70,000 computational cycles, this step is performed to harmonize the water molecules with the protein structure. The third step is the final molecular dynamics step for 120 ns corresponding to 60,000,000 computational cycles (38). The time step was set at 2 fs and the simulations were performed under NpT conditions at both constant temperature (310 K) and pressure (1 atm), using Langevin dynamics with a damping constant of 1 ps-1. The rmsd, total energy per frame and hydrogen bonding were calculated directly by existing commands in the analysis menu of the VMD software while the centre of the mass radius of gyration was calculated by VMD using the Tkconsole tool and based on previously prepared scripts.

Molecular mechanics Poisson-Boltzmann surface area method (MM/PBSA).
In this work, we were interested in determining the relative binding energy in order to be able to judge the efficiency of both ligands (azithromycin, dirithromycin) to block the targeted protein. This calculation was performed using an end-point approach called the molecular mechanics Poisson-Boltzmann surface method (MM/PBSA) via the opensource program CaFE (Calculation of Free Energy) (39).
An MM-PBSA script was prepared with the following set of parameters, an outer dielectric constant was set to 80.0, the inner dielectric constant to 1.0 and the grid spacing of 0.05 Å was employed, while for the SA calculation, the surface tension value was fixed at 0.00542 with a surface offset of 0.92. Finally, the binding energy was averaged over a set of conformations (100 snapshots) as follows (40).

Results and Discussion Molecular docking study
This work was carried out to predict the potential activity of the antibiotic dirithromycin aiming to assess its therapeutic potential in relation to COVID-19. The results were compared to those of the wellknown azithromycin which is a molecule highly described also for its anti-viral potential All target molecules have been subjected to three-dimensional modeling and molecular docking with three families of "protein" receptors involved in the mechanism of infection of the SARS-COVID-19 virus including (i) three cell host binding proteins (Furin, ACE2 and TMPRSS-2); (ii) eight SARS-CoV (Receptor binding domain of Spike protein, NSP9 RNA binding protein, NSP3 ADP ribose phosphatase, RNA binding domain, NSP15 endoribonuclease, nucleocapsid protein Nterminal RNA binding domain, papain-like protease and 3CL-protease); and (iii) four pro-inflammatory mediators (COX2, PLA2, NIK and IRAK-4). As shown from the results depicted in Supplementary materials (S1 and S2), the two tested ligands displayed similar affinity towards the different selected target proteins. However, despite the similarity of the affinity between the two ligands with respect to the different target receptors, there is a difference in the level of interaction energy. In fact, dirithromycin possesses the lowest binding energy towards all receptors compared to azithromycin, which gives it a greater potential for inhibition. Interestingly, the docking results allowed us to classify the obtained complexes into five different families with different energy levels.
The first class exhibited the lowest binding energy between the ligands and both furin and TMPRSS-2 receptors. Dirithromycin binds better to furin (-9.9 Kcal/mol) than azithromycin (-9.4 Kcal/mol). Regarding the mode of their interaction with the catalytic site, dirithromycin can interact and block two amino acids of the furin catalytic triad, namely Asp154 and Ser368 by van der Walls interactions, and can also interact with two amino acids of the same triad by hydrogen bonds, namely His194 and Ser368, while azithromycin interacts only with van der Walls bonds between the catalytic triad and only Asp154 and His194 residues and likewise with the metalbinding site Asp258 via van der Walls interactions. Alike, there is a difference in the binding energy of dirithromycin-TMPRSS-2 complex (-10.7 Kcal/mol) and azithromycin-TMPRSS-2 (-9.0 Kcal/mol) whilst no ligand-amino acid interactions in the active site have been established. Consequently, dirithromycin may be potent furin and in fewer degree TMPRSSinhibitors.
The second family implies the two interaction complexes dirithromycin and azithromycin with the RNA binding domain (PDB ID: 6VYO). In this case, dirithromycin showed a slightly high affinity for the receptor with interaction energy of -9. 7 Kcal/mol at the RNA binding domain compared to that of azithromycin, -9.5 Kcal/mol. The SARS-CoV-RNA binding domain protein is one of the most attractive targets for the development of new drugs against CoVs because of its important role in the replication and transcription of the virus. This structure includes the RNA binding domain (amino acids 50-173) which can be a promising target to control or stop the spread of the virus. Based up on our docking study, we noticed that dirithromycin strongly interacts with this protein by binding in the internal space resulting from the tetrameric formation of this protein. This ligand allowed us to obtain six solutions of dirithromycin-RNA complex (binding 6VYO) with binding energy values from -9. 7 Kcal/mol to -9.0 Kcal/mol, while the results obtained with azithromycin have shown ligand-receptor complexes with energies in the range -9. 5 Kcal/mol to -8. 6 Kcal/mol.
The last family includes the RBD of the protein Spike, the NSP9, NSP3 ADP ribose phosphatase, NSP15 Endoribonuclease, Nucleocapsid N-terminal protein RNA binding domain, NIK and PLA2 where the interaction energies of the two ligands are very high and do not predict the binding of the two ligands to the proteins described above. However, it is important to note that dirithromycin retains the highest affinity with respect to the various receptors compared to azithromycin, except with NSP3 ADP ribose phosphatase complexes where the azithromycin-NSP3 complex has present higher binding energy (-9.3 Kcal/mol) than dirithromycin-NSP3 complex (-8. 6 Kcal/mol), whereas dirithromycin-NSP15 Endoribonuclease and azithromycin-NSP15 Endoribonuclease exhibited equipotent binding energies (-8.4Kcal/mol).
The COX2-Ligand and ACE2-Ligand interactions are part of the fourth family of complexes ( Table 2). Analysis of the ACE2-Ligand interaction complex shows that dirithromycin binds to two amino acids at the metal-binding site, His378 by Pi-sigma bond and Glu402 by van der Walls bond, respectively, while azithromycin binds to the same amino acids only with van der walls bonds. The binding energy levels are also in favor of dirithromycin with a difference equal to 1.3 Kcal/mol (dirithromycin-ACE2 -10.6 Kcal/mol vs. azithromycin-ACE2 -9. 3 Kcal/mol Additionally, previous results have demonstrated that teicoplanin, an anti-Gram-positive glycopeptide drug, could be used as an alternative medication against the novel SARS/COVID-19 (7). This antibiotic controls also various viruses including some coronavirus such as Middle East respiratory syndrome coronavirus (MERS-CoV) and SARS-CoV (15)(16)(17). It has been demonstrated that this drug inhibits the cleavage of the spike protein at low pH in SARS-Cov-2 (15). Bafilomycin A1 is another macrolide considered a promising anti-COVID-19 drug. It is an endo/lysosomal V-ATPase inhibitor that can stop SARS-CoV and SARS-CoV-2 virus by inhibiting the function of the ACE2 receptor (19).
Using the VMD software, the calculation of the RMSD and the RMSF was carried out. RMSF was measured relative to the C atom of each amino acid for the protein in complexes with azithromycin or dirithromycin as well as the free protein. The superimposed RMSF plot (Figure 4) showed different ranges of variation where it is noted that the most significant fluctuations were observed at the level of the curve representing the free protein with a residual fluctuation which varies between 16.8 Å and 19.9Å. At the level of the furin-dirithromycin complex, a decrease in fluctuation was observed to stabilize in a range of variation between 13.2 Å and 17.5Å while the minimum fluctuation was observed from the trace which present the protein in complex with azithromycin and fluctuation values of 10.7Å -13.8Å have been observed.
The calculation of the root mean square deviation (RMSD) provides information on the stability of the complexes by comparing it to the RMSD of the protein alone ( Figure 5). The superposition of the RMSD traces allowed us to determine the first part during the first 48 ns where the protein alone, in complex with azithromycin and the protein in complex with dirithromycin, exhibit the same degree of stability which is about 1.5 Å. Beyond 48ns, the stability of the free protein, as well as the protein in complex with azithromycin, decreases throughout the simulation to stabilize at values of 2.5 Å and 2 Å respectively while the protein in complex with the dirithromycin exhibit a stable plateau throughout the simulation with RMSD values under than 1.7Å. These results lead to the conclusion that both ligands provide stability by binding to the protein, with advantageous stability for dirithromycin compared to azithromycin targeting the furin protein.
Calculating the radius of gyration (Rg) allows us to extract additional information regarding the compactness of a protein as well as its folding properties. Low Rg values are the result of a tight three-dimensional structure while high Rg values prove the conformational instability of a protein (29). By comparing the different Rg plots obtained in Figure 6, during the first 60ns of simulation no difference was observed, the three Rg plots of the free protein, furin-azithromycin and furin-dirithromycin were stable at values in order of 22.2Å, then a slight increase in the Rg values was observed at the level of the graph representing the furin-azithromycin complex to stabilize at values of 22.5Å. This increase was not observed in the furin-dirithromycin complex which exhibits stable values and does not exceed 22.25Å. Moreover, even if the Rg results do not show a large difference, however, they are in favor of dirithromycin which offers greater compactness compared to azithromycin which decreases the compactness of the protein even in comparison with the free protein. Unlike this, it is important to monitor the number of hydrogen bonds formed throughout the simulation which play a very important role in the formation and stability of the ligand-receptor complex. It should be noted that the more hydrogen bonds the protein has, the more stable its conformational structure is during the simulation (44). In our work, the number of total hydrogen bonds formed at the protein level was calculated per frame throughout the simulation (Figure 7). A first analysis was carried out to compare the number of hydrogen bonds formed at the level of the complexes compared to the free protein ( Figure 7A). Figure 6. Analysis of Rg plots of furin-free protein (black) and complexed protein with azithromycin (red) and dirithromycin (bleu) at 120 ns MD simulations. Figure 7. Analysis of the total number of H-bond count throughout the simulation of furin free protein (black), a complexed protein with azithromycin (red) and dirithromycin (bleu) at 120 ns MD simulations (A), Difference in H-bond number between furin-azithromycin complex-free protein (red) and furin-dirithromycin complex-free protein (bleu) at 120 ns MD simulations (B). The difference in the number of hydrogen bonds formed between dirithromycin and azithromycin throughout the simulation is shown in C.
An increase in the numbers of hydrogen bonds compared to the free protein following the binding of the two ligands azithromycin and dirithromycin was noticed (Figure 7 B and C). The average of the hydrogen bonds per frame at the level of the free protein is 120 H-bonds per frame followed by an average of 122 H-bonds per frame at the level of the protein in complex with azithromycin and a higher average at the level of the furin-dirithromycin complex which has an average of 129 H-bonds per frame. Comparing azithromycin and dirithromycin for hydrogen bonds formation, we notice a difference of 7 H-bonds per frame which is well visualized in Figure  8. This leads us to better understand the effectiveness of dirithromycin to block the target protein compared to azithromycin.
The analysis of the energetic behavior was carried out for the two complexes using the software VMD where we started to follow the electrostatic energy and van der Walls to have total energy calculated per frame. This was superimposed on the plot of center mass (com) calculated as described previously to determine the distance separating the center of mass of the ligand from the center of mass of the protein during the simulation to judge the behavior of the ligand throughout the simulation. In Tables 3 and 4, we notice that at the level of the furin-azithromycin complex, the total energy varies at the start of the simulation and has values between -40 kcal/mol and 0 kcal/mol to stabilize at the 48 ns with values that vary between -20 and -10 kcal/mol for then passed to a plateau of 0 kcal/mol observed at the level of the last 12 ns of the simulation.
The same behavior was observed at the level of the com's ( Figure 8) where we observe a detachment of the ligand at the level of the first part of the simulation with values of com's in the order of 50 Å to then stabilize at values in the range of 25 Å and during the last 12 ns of the simulation, the values increase to reach 60 Å.   -9783.49 -1939-9783.49 - .91 -212.97 108.76 -11723.40 -104.20 -9996.46 -1831-9783.49 - .14 -11827.61 Receptor -9777.22 -1918-9783.49 - .57 -212.96 107.45 -11695.79 -105.50 -9990.18 -1811 ε Elec: electrostatic energy change; ε Vdw: van der Waals energy change; ε PB: Poisson-Boltzmann; ε SA: Surface Area; ε Gas: gas-phase molecular mechanical energy change; ε Sol: the solvation free energy change; ε Pol/ ε NPol: polar and non-polar contributions to the solvation free energies.

Figure 8.
Graphical superposition of the trace presenting the distance between both centers of mass furin and azithromycin with the variation of the total energy during 120 ns (A) and the distance between both centers of mass furin and dirithromycin with the variation of the total energy during 120 ns (B).
From these analyzes, it was possible to identify the source of the instability energy of azithromycin which is due to the detachment of the ligand. While at the level of dirithromycin and under the same simulation conditions, no detachment was observed. The total energy values vary between -80 kcal/mol and -40 kcal/mol to stabilize at the level of the last 24 ns at -20 kcal/mol. Similarly, the absence of aberrations, which are signs of detachment, can be seen in the com's plot, with values varying between 18-28 Å.
To decide on the affinity of the two ligands and their efficiency in blocking the furin protein, we have calculated the relative binding energy MM-PBSA which is one of the most reliable and widely used approaches and which combines the models of molecular mechanics and continuous solvent to calculate G-bind of small molecules (44)(45)(46). By comparing the total relative binding energy of the furin-azithromycin complex which is -18.2793 kcal/mol compared to the relative binding energy presented by the furin-dirithromycin complex which is -51.6722 kcal/mol we can conclude on the effectiveness of dirithromycin in blocking the furin with energy 2.8 times higher than that observed at the level of azithromycin under the same working conditions.

Conclusions
In this study, we have explored for the first time the re-purposing of existing drugs dirithromycin as anti-SARS-CoV-2. Our molecular docking study suggests that dirithromycin appears to be the most powerful SARS-CoV-2 inhibitor than azithromycin exhibiting good binding affinity with cell host binding proteins, especially with furin and TMPRSS-2 and displaying potential anti-inflammatory properties. MD simulations explain and confirm why dirithromycin is more effective than azithromycin in treating COVID-19, especially in blocking furin with an energy level 2.8-fold higher. Furthermore, dirithromycin should be considered as an alternative for the treatment of the SARS-CoV 2 virus by evaluating its effectiveness using both in vitro and in vivo experiments.

Funding
This research has been funded by the Scientific Research Deanship at the University of Ha'il-Saudi Arabia through project number RG-20141.