key: cord-0790959-e6sgloe8
authors: Al-Bustany, Hazem Abbas; Ercan, Selami; Ince, Ebru; Pirinccioglu, Necmettin
title: Investigation of angucycline compounds as potential drug candidates against SARS Cov-2 main protease using docking and molecular dynamic approaches
date: 2021-04-10
journal: Mol Divers
DOI: 10.1007/s11030-021-10219-1
sha: 7607ec7d867c143be04246d1e9ac67daa73779d8
doc_id: 790959
cord_uid: e6sgloe8

ABSTRACT: The emerged Coronavirus disease (COVID-19) causes severe or even fatal respiratory tract infection, and to date there is no FDA-approved therapeutics or effective treatment available to effectively combat this viral infection. This urgent situation is an attractive research area in the field of drug design and development. One of the most important targets of SARS-coronavirus-2 (SARS Cov-2) is the main protease (3CLpro). Actinomycetes are important resources for drug discovery. The angucylines that are mainly produced by Streptomyces genus of actinomycetes exhibit a broad range of biological activities such as anticancer, antibacterial and antiviral. This study aims to investigate the binding affinity and molecular interactions of 157 available angucycline compounds with 3CLpro using docking and molecular dynamics simulations. MM-PBSA calculations showed that moromycin A has a better binding energy (− 30.42 kcal mol(−1)) compared with other ligands (in a range of − 18.66 to − 22.89 kcal mol(−1)) including saquayamycin K4 (− 21.27 kcal mol(−1)) except the co-crystallized ligand N3. However, in vitro and in vivo studies are essential to assess the effectiveness of angucycline compounds against coronavirus. GRAPHIC ABSTRACT: [Image: see text] SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s11030-021-10219-1.

Wuhan Novel Coronavirus (2019-nCoV), which is known as Coronavirus disease 2019 (COVID-19) outbreak, has become worldwide pandemic and caused deaths of numerous infected people [1] . So there is an immediate and a large-scale research campaign to design and discover definitive therapeutic agents against this fatal virus. To date, there is no effective antiviral therapy [2] , but few vaccines have recently been approved to combat the COVID-19 infection. Most of the immediate efforts focused on clinically approved known drug evaluation and virtual screening of available molecules from chemical libraries [3] .

The main protease (3CLpro, also called M pro ), which is involved in processing the polyproteins that are translated from the viral RNA [4] , is one of the most sought drug targets among coronaviruses [5] . There are no human proteases with similar cleavage specificity. Therefore, inhibitors for this protein are expected not to be competitive for the host.

There have been recent studies on developing inhibitors targeting 3CLpro based on structure-based ligand design and molecular modelling and even in vitro studies from available libraries [6] [7] [8] . This study aims to investigate the potentials of angucycline compounds as drug candidates against coronavirus. The angucyclines are the largest group of polyketides produced in streptomycetes via complex enzymatic biosynthetic pathways [9] and show a broad range of biological activities, such as antibacterial [10, 11] , antitumor [12] , antiviral [11, [13] [14] [15] and enzyme inhibitors [11] . The binding affinity of a total of 157 angucycline compounds against 3CLpro was screened by molecular docking. Molecular dynamics (MD) simulations were performed for the complexes of 3CLpro with saquayamycin K4 and moromycin A, which produced the best docking scores and also for some other crystal 3CLpro complexes, namely N3 [7] , 11r, 13a, 13b and 14b [6] . The binding energies for the complexes were calculated by MM-PBSA.

The crystal structure (PDB ID: 6LU7) for 3CLpro is obtained from Protein DataBank (PDB) (http:// www. rcsb. org) [6] . All the angucyclines were previously prepared and ready for docking study [16] . They were drawn by ChemDraw and converted to 3D structures and minimized by Discovery Studio Visualizer [17] . The ligands and the protein were prepared by MGL Tools program [18] , and they were saved in PDBQT file format. The binding site of 3CLpro enzyme was selected for grid calculations by taking the advantage of the co-crystallized N3 ligand contained in the crystal structure of 6LU7.pdb file. A grid box with the dimensions of 18 × 24 × 20 Å was placed to − 9.66, 11.50 and 68.81 points through x, y and z coordinates, respectively. Each running gave a total of 9 docking poses for each ligand. Molecular docking was performed using AutoDock Vina 1.1.2 [19] , running on Windows 8.1 operating system with 4 CPUs and 12 GB of RAM that generated 9 poses for each docked ligand.

The ADME properties of saquayamycin K4 and moromycin A were defined by using SwissADME Web site [20] .

All MD simulations were conducted using the Assisted Model Building with Energy Refinement (AMBER, v.19.0) software package [19] running on TRUBA clusters (TUBI-TAK). AM1-Bcc (the Austin model with bond and charge correction) atomic partial charges for the ligands were calculated by the antechamber module of the AMBER (v.19.) package [21] . Xleap as implemented in AMBER was employed to prepare parameter/topology and coordinate files and to solvate and neutralize the system for the MD simulations. The complexes were solvated in an octagonal periodic box containing explicit TIP3P [22] with dimensions that its boundaries extend at least 10 Å from any solute atom. This is a very popular solvation model for MD simulations because of its simplicity and computational efficiency. The ff19SB force field was used for the protein, while the General Amber Force Field (GAFF) [23] was employed for the ligands because it can handle small organic molecules. Three-dimensional and surface structures were displayed using the Chimera (UCSF) software package [24] .

The complex including water and sodium ions was minimized in two steps. In the first step, the enzyme and the ligand were kept fixed; only the water molecules were allowed to move. In the second step, all of the atoms were allowed to move. In both steps, energy minimization was performed in 2000 steps: 1000 steps with the steepest descent method and 1000 steps with the conjugate gradient method. Heating was performed for 200 ps, where the protein-ligand complex was restrained with a force constant of 50 kcal mol −1 Å −2 . Equilibration was performed for 1 ns in the isothermal-isobaric (NPT) ensemble, restraining the protein-ligand complex by 1 kcal mol −1 Å −2 . Final simulations (the production phase) were performed for 50 ns in the NPT ensemble at a temperature of 300 K and a pressure of 1 atm. The step size was 2 fs during the entire simulation. A Langevin thermostat and barostat were used to couple the temperature and pressure. The SHAKE algorithm was applied to constrain all bonds containing hydrogen atoms [25] . The nonbonded cutoff was kept at 10 Å, and longrange electrostatic interactions were treated using the particle mesh Ewald (PME) [26] method with a fast Fourier transform grid with a spacing of approximately 0.1 nm. Trajectory snapshots were taken 1-ps intervals, and these were later used in the analysis.

The root-mean-square deviation (RMSD) and rootmean-square fluctuation (RMSF) were analyzed by cpptraj as implemented in AMBER from MD simulations, and the results were plotted using MS Excel.

The Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) or Molecular Mechanics Generalized Born Surface Area (MM-GBSA) module of AMBER (v.19.0) was applied to compute the binding free energy (ΔG bind ) of each complex [27, 28] . The application of the method is detailed in our previous works [29, 30] . For each complex, 2500 frames were extracted from a total of 25,000 frames calculated over a period of 50 ns.

The 2D structures of the compounds saquayamycin K4, moromycin, N3, 11r, 13a, 13b and 14b are shown in Fig. 1 . The structures of all the angucyclines are provided in Supplementary Information (Table S1 ). The active site of 3CLpro consists of a large pocket with deeper three wells ( Fig. 2) , involving a catalytic dyad composed of CYS145 and HIS41 and four subsites (S1, S2, S3 and S4) [31] [32] [33] [34] [35] The S1 subsite involves PHE140, LEU141, ASN142, HIS163, GLU166 and HIS172 residues. This subsite is further divided which includes THR25, THR26 and LEU27. The S2 subsite has a hydrophobic character, which is comprised of HIS41, MET49, TYR54, MET165 and ASP187. On the other hand, the S3 subsite is solvent-exposed, and therefore, it can accommodate a wider range of functional groups [6, 8] . The S4 binding subsite is composed of MET165, LEU167, PHE185, GLN189 and GLN192 residues [6] .

To validate the application of AutoDock Vina for the docking of the ligands into the protein, N3 ligand was extracted and re-docked into the active site of the enzyme and the generated docked poses were compared with the co-crystallized pose with a RMSD of 1.711 Å. Docking of N3 by Dock 6 [36] also produced a similar conformational pattern with the co-crystallized ligand with a RMSD of 1.390 Å (Fig. 3) .

All the compounds were well docked into this pocket using AutoDock Vina software (Figs. 4, 5, 6, 7, 8, 9, 10) . Docking analyses indicated that 12 candidates have docking scores lower (the lower the score the better the binding affinity) than − 9.5 kcal/mol as summarized in Table 1 . They all have better docking scores compared to the current inhibitors, N3 [7] , 11r, 13a, 13b and 14b [6] . The top four candidates to bind with 3CLpro among the angucycline compounds were saquayamycin K4, moromycin A, saquayamycin B and saquayamycin B-2. The docking scores of the rest of angucyclines are demonstrated in Supplementary Information (Table S2 ). They all bind to the protein via hydrogen bonds and hydrophobic interactions.

MD simulation was performed for the complexes of the best scoring ligands, namely saquayamycin K4, moromycin A and the other ligands (N3, 11r, 13a, 13b and 14a); each complex was simulated for a period of 50 ns. The MD simulations of these two angucyclines are expected to reflect the mode of the action of the rest of angucyclines. The rootmean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) obtained from MD trajectories for the complexes of 3CLpro with all the ligands are presented in Figs. 4 and 5. The lowest energy structures for the complexes are demonstrated in Figs. 6, 7, 8, 9, 10, 11, 12 as we will see in the discussion section. The binding energies calculated from MD trajectories by MM-PBSA are summarized in Table 1 .

The COVID-19 has caused one of the deadliest outbreaks in the present time. Prevention of the newly emerging COVID-19 infection is very challenging. The possibility of in silico methods can be utilized to find solutions with fewer trials and errors and thus saving both time and costs of the research. The docking study provides a time saving, shorter and cheaper approach for drug development.

RMSD and RMSF obtained from MD trajectories for the complexes of 3CLpro with all the ligands (Figs. 4 and 5) demonstrate that all the complexes have reached a conformationally stable structure and the residues involved in the interaction with the ligands, and ligands themselves are in compact status and they do not undergo significant dynamical changes during MD simulations except saquayamycin K4.

The comparison of the superimposed structures of the docked N3 by AutoDock Vina with the original X-ray coordinates indicates that the application of AutoDock for the docking of the ligands into the protein is reliable (Fig. 3) . The N3 inhibitor has relatively a lower docking score (− 7.8 kcal/mol) compared with almost of 2/3 of the angucycline compounds (Table S1 ). The surface filling model shows that this ligand fully occupies the cavity in the active site (Fig. 6) . It also involves hydrogen bonds and hydrophobic interactions to bind with 3CLpro (Fig. 6) .

After the MD calculation, the ligand has not undergone significant conformational changes within the active site compared to saquayamycin K4 and moromycin A. Although N3 is fitted well to the active site, it does not have so specific interactions compared to Moromycin A. The main noticeable interactions are that 1,2-isoxazole ring is held within a cavity composed of PRO168, THR190 and ALA191 residues via CH-pi interactions, and the benzyl ring is sandwiched between HIS41 and MET49 via pi-pi and CH-pi interactions (Fig. 6) . The CYS44 residue has also an interaction with the benzyl ring via pi-S interaction. The lactam ring interacts with SER144 and CYS145 via hydrogen bonds, and more markedly, GLN189 interacts with alanine and GLU166 with leucine next to alanine in N3 via hydrogen bonding (Fig. 6) . Docking results showed that saquayamycin K4 has the highest docking score against 3CLpro among angucycline compounds (Table 1 ) with a score of − 11.0 kcal/mol. The ligand is well located in the active site, filling all the area in the pocket (Fig. 7) . The binding mode of this compound involves hydrogen bonds and hydrophobic interactions. The complex of the compound with the protein includes hydrogen bonds with MET49, ASN142, GLY143, CYS145, HIS164, GLU166 and THR190 and hydrophobic interactions with LEU27, MET49 and CYS145. The binding also involves an NH-pi stacking with GLN189. HIS164, together with HIS41 and CYS145 residues, is located in the active site of the protein, and they are directly involved in the catalytic cleavage process. It is also reported that GLN189 is also involved in the catalytic process [8] . The indole ring of saquayamycin K4 has a crucial role in the binding, occupying one of three wells in the active site (Fig. 7) , thus contributing to the binding energy and blocking the role of GLN189 in the catalytic process. On the other hand, the methyl group on the H ring of Saquayamycin K4 is also cradled in one of the other wells and interacts with LEU27 via hydrophobic interactions. After MD simulations, the ligand underwent a significant conformational change within the active site (Fig. 7) . The binding energy obtained from MD simulations demonstrates that this compound has a comparable binding energy to the rest of the ligands except N3 (Table 1 ). The lowest energy structure for the complex of saquayamycin K4 with 3CLpro from MD calculations indicates that the H and I rings of the ligand lay out of the active pocket; however, the rest covers the major surface area of the active site, in a manner of "stone cave door" (Fig. 7) . In this way, the residues involved in the catalysis are totally expected to be blocked from being accessed by the native substrate. The A, B, C and indole rings are held in one of the active site pockets, mainly made of hydrophobic residues. These two parts interact with LEU167, ALA191, PRO168 and the hydrophobic portion (-CH 2 CH 2 -) of GLU166 and GLN189 on their side chains where GLN189 is thought to hold a meaningful role in the catalytic activity of the protein [8] . This residue as well as SER46 and GLU47 is involved in a hydrogen bond network with the ligand. GLU47 has twopoint interactions with the hydrogen donors attached to two Fig. 2 Surface representation of the active site of 3CLpro highlighting its sub-pockets. Yellow and green represent the S1 subsite, blue represents the S2 subsite and red represents the S4 subsite hydroxyls linked to the bridged carbons between the F and G rings, marked as H12 and H13 (Fig. 1) , while SER46 has similarly two-point contacts with the phenolic oxygen (O20) and the carbonyl (O21) on the D ring. On the other hand, the amide hydrogens of GLN189 interact with the oxygen acceptors marked as O23, O26 and O38. So the hydroxyl group bearing H13 and O26 is chelated between GLU47 and GLN189. It is possible that this combination of the hydrogen bond network involving SER46, GLU47 and GLN189 contributes to a stronger binding of this ligand to the protein.

Moromycin A has a docking score of − 10.8 kcal/mol ranking the second with a similar binding mode of saquayamycin K4 (Fig. 8) . It forms hydrogen bonds with the residues GLY143, CYS145, GLN189, and THR190, hydrophobic interactions with the residues LEU27, CYS145, PRO168 and ALA191 and alkyl-pi interactions with the residues HIS41, HIS163, MET165 and PRO168. In this complex, moromycin A also interacts with the main residues involved in the catalytic process. The binding energy obtained from MD simulations demonstrates that this angucycline has relatively comparable binding energy compared to the cocrystallized N3 ligand (Table 1 ) but higher binding energies compared to the other ligands (11r, 13a, 13b and 14b11r,  13a, 13b and 14b) . The lowest energy conformer obtained from MD simulations for the complex of Moromycin A with 3CLpro indicates that the ligand resides in the active pocket composed of all the four sub-pockets (Fig. 8 ). The C, E and F rings of the ligand are settled within the hydrophobic pocket of the active site. C ring interacts with ALA191 and ALA193 residues via van der Waals interactions, whereas E and F rings interact with LEU167, PRO168, MET49 and MET165 via CH-pi interactions. The G ring also interacts with MET49 and MET165 via van der Waals interactions. The carbonyl groups on the E and G ring form hydrogen bonds with GLN189 and 192 (Fig. 5) . On the other hand, the methyl and CH 2 groups in the H ring have CH-pi interactions with the crucial residue HIS41, and the I ring is held within a pocket consisted of THR25, THR26, LEU27, VAL42 and CYS145 via both hydrophobic and hydrogen bonding (Fig. 8) . In a similar pattern to N3 inhibitor, 11r, 13a, 13b and 14b filled the active pocket in 3CLpro (Figs. 9, 10, 11, 12) with quite close docking scores. Following the MD simulations, all the ligands kept their conformation in the active sites with only slight changes compared to saquayamycin K4 (Figs. 9, 10, 11, 12 ). 11r has hydrogen bonds with CYS145 and HIS163 via the lactam ring and 3-phenylpropanoyl residue with GLN189. The cyclohexyl ring is located in one of the cavities in the pocket and N-benzyl group interacts with MET49 and CYS44 via CH 2 -pi and S-pi interactions, respectively. Yet it is interesting to see that one molecule of water has penetrated into the cavity residing between HIS41, CYS145, and the carbonyl group of the lactam ring (Fig. 9) . The cyclohexyl group and the lactam ring in 13a are located in the active site in a similar manner to 11r where the tert-butoxy carbonyl residue is located within a hydrophobic pocket consisted of MET49, THR24, THR25 and THR45. Interestingly, one water molecule has penetrated in the same cavity as in 11r, resided between HIS41 and ASP187 where the cyclohexyl ring is sandwiched between MET49 and 165 (Fig. 10) via hydrophobic interactions. This ring has also a contact with HIS41 via CH-pi interactions. Moreover, the cyclopropyl ring interacts with PRO168 and ALA191 via hydrophobic interactions. Lastly, the critical residue GLN189 has a hydrogen bond with NH attached to the cyclopropyl. In the complex of 13b, the location of the lactam ring is switched with N-benzyl group compared to those in 11r and 13a and the cyclopropyl ring is buried in the location where the cyclohexyl as in 11r and 13a (Fig. 11) . The significant interactions are the hydrogen bonds of lactam ring with ASN119 and CYS145 residues. The tert-butyl group interacts with PRO168 via hydrophobic interactions and benzyl group interacts with ASN142 via NH-pi interactions. Furthermore, the cyclopropyl ring is sandwiched between HIS41 and MET165 via CH-pi and hydrophobic interactions, respectively. In this case, two molecules of waters were found in the active site one interacts with the carbonyl group on 3-amino-1H-prydine-2-one ring, and the other interacts with the carbonyl group close to the cyclopropyl ring. The lowest energy structure of 13b obtained from MD simulations indicates that the orientation of the groups in the active sites resembles that of 14b where 3-amino-1H-prydine-2-one interacts with GLU166 and HIS41 interacts with 3-cyclopropanoyl carbonyl (Fig. 11) . Finally, as to the evaluation of the complex of 14b with 3CLpro, it is evident that 14b has a similar binding mode to that of 13b with regard to the orientation of the groups in the active site (Fig. 12) . The lactam ring has hydrophobic interactions with THR25, LEU27 and VAL42, while HIS41 has hydrogen bonds with 3-cyclopropylpropanoyl carbonyl, and GLU166 and GLN189 have hydrogen bonds with the amino and carbonyl groups in 3-amino-1H-prydine-2-one ring, respectively. On the other hand, the interactions of cyclopropyl ring with MET165 via hydrophobic interactions, 3-amino-1H-prydine-2-one ring with MET49 via CH-pi interactions and the benzyl group with ASN142 via NH-pi interactions also contribute to the binding of this ligand to the protein. The above findings indicate that moromycin A and saquayamycin K4 showed high affinity and good binding interactions. These compounds have potential inhibitory effects for the catalytic dyad Cys145 and His41 as well as the other amino acids involved in the binding pocket. So their ability for the interaction with 3CLpro suggests further gains in the inhibition of the virus activity. Moreover, the interaction patterns of these compounds possess advantages compared to the known FDA drugs. This may make them potential for the treatment of COVID-19.

The drug-likeness properties of moromycin A and saquayamycin K4 were calculated by SwissADME. The results showed that both ligands have two violations according to Lipinski's rule of five due to their masses and containing more than 10 hydrogen bond acceptors. While both ligands seem to demonstrate P-glycoprotein substrate properties, saquayamycin K4 seems also to inhibit Cytochrome P450 2C9 enzyme. Fig. 9 (a) 11r docked to 3CLpro, (b) 2D structure of 11r docked to 3CLpro and (c) the binding mode of the complex of 11r with 3CLpro in the lowest energy structure obtained from MD simulations Fig. 10 (a) 13a docked to 3CLpro, (b) 2D structure of 13a docked to 3CLpro and (c) the binding mode of the complex of 13a with 3Clpro in the lowest energy structure obtained from MD simulations

The protein-ligand complex interactions play a significant role in structure-based drug design. Docking analyses showed that the angucycline compounds are interacted with 3CLpro active site by hydrogen bonds and hydrophobic interactions, almost 2/3 of them with relatively higher scores compared to some reported inhibitors. MD simulations and MM-PBSA results for the complex of saquayamycin K4, moromycin A and inhibitors indicated that moromycin A has a similar binding energy compared to that of the co-crystallized N3 inhibitor. Therefore, this compound and the others with the similar structural properties are considered promising lead candidates against COVID-19. The structure of moromycin A may be also used as a template for the design of new drug candidates against 3CLpro. However, in vitro and in vivo works are required to achieve a final conclusion about their efficiency as anti-coronavirus candidates. 

The Essential Facts of Wuhan Novel Coronavirus Outbreak in China and Epitope-based Vaccine Designing against COVID-19

Rapid identification of potential inhibitors of SARS-CoV-2 main protease by deep docking of 1.3 billion compounds

Ivanenkov Y (2020) Potential COVID-2019 3C-like protease inhibitors designed using generative deep learning approaches

From SARS to MERS: crystallographic studies on coronaviral proteases enable antiviral drug design

Coronavirus main proteinase (3CLpro) structure: basis for design of anti-SARS drugs

Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors

Identification of potential binders of the main protease 3CLpro of the COVID-19 via structure-based ligand design and molecular modeling

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α-ketoamide inhibitors

Crystal structures of two aromatic hydroxylases involved in the early tailoring steps of angucycline biosynthesis

Angucyclines: total syntheses, new structures, and biosynthetic studies of an emerging new class of antibiotics

Angucycline group antibiotics

Mechanisms underlying the anticancer activities of the angucycline landomycin E

Secondary metabolites by chemical screening. Part 19. SM 196 A and B, novel biologically active angucyclinones from Streptomyces sp

New antiviral antibiotics, cycloviracins B1 and B2. I. Production, isolation, physicochemical properties and biological activity

New antiviral antibiotics, cycloviracins B1 and B2. II Structure determination

Screening of angucycline antibiotics as potential drug candidates against MRSA by docking analysis

Discovery Visualizer Studio

AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility

SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method

Comparison of simple potential functions for simulating liquid water

Development and testing of a general amber force field

UCSF Chimera-a visualization system for exploratory research and analysis

A second generation force field for the simulation of proteins, nucleic acids, and organic molecules

Particle mesh Ewald: an N⋅log(N) method for Ewald sums in large systems

Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models

Continuum solvent studies of the stability of DNA, RNA, and PHOSPHORAMIDATE−DNA helices

Experimental and theoretical study of the mechanism of hydrolysis of substituted phenyl hexanoates catalysed by globin in the presence of surfactant

Computational design of a fulllength model of HIV-1 integrase: modeling of new inhibitors and comparison of their calculated binding energies with those previously studied

Structural basis for the inhibition of SARS-CoV-2 main protease by antineoplastic drug carmofur

Tackling COVID-19: identification of potential main protease inhibitors via structural analysis, virtual screening, molecular docking and MM-PBSA calculations

Identification of potential molecules against COVID-19 main protease through structure-guided virtual screening approach

Fast identification of possible drug treatment of coronavirus disease-19 (COVID-19) through computational drug repurposing study

Identification of potential Mpro inhibitors for the treatment of COVID-19 by using systematic virtual screening approach

DOCK 6: impact of new features and current docking performance

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgements The authors thank The Scientific and Technological Research Council of Turkey (TUBITAK) for providing TRGRID facilities and DA Case (University of California, San Francisco) for a waiver license for AMBER (v.19.0).

The authors declare that they have no conflict of interest.

Hazem Abbas Al-Bustany 1 · Selami Ercan 2 · Ebru Ince 3 · Necmettin Pirinccioglu 4 1