key: cord-0996824-64yueq69
authors: Almeida, Joyce S. F. D.; Botelho, Fernanda D.; de Souza, Felipe R.; dos Santos, Marcelo C.; Goncalves, Arlan da Silva; Rodrigues, Rodrigo L. B.; Cardozo, Monique; Kitagawa, Daniel A. S.; Vieira, Leandro A.; Silva, Raphael S. F.; Cavalcante, Samir F. A.; Bastos, Leonardo da C.; Nogueira, Mariana de O. T.; de Santana, Priscila I. R.; Brum, Juliana de O. C.; Nepovimova, Eugenie; Kuca, Kamil; LaPlante, Steven R.; Galante, Erick B. F.; Franca, Tanos C. C.
title: Searching for potential drugs against SARS-CoV-2 through virtual screening on several molecular targets
date: 2021-01-08
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2020.1869096
sha: c6ddba65256f3a9f455fa7372c97ae4e33434888
doc_id: 996824
cord_uid: 64yueq69

The acute respiratory syndrome caused by the SARS-CoV-2, known as COVID-19, has been ruthlessly tormenting the world population for more than six months. However, so far no effective drug or vaccine against this plague have emerged yet, despite the huge effort in course by researchers and pharmaceutical companies worldwide. Willing to contribute with this fight to defeat COVID-19, we performed a virtual screening study on a library containing Food and Drug Administration (FDA) approved drugs, in a search for molecules capable of hitting three main molecular targets of SARS-CoV-2 currently available in the Protein Data Bank (PDB). Our results were refined with further molecular dynamics (MD) simulations and MM-PBSA calculations and pointed to 7 multi-target hits which we propose here for experimental evaluation and repurposing as potential drugs against COVID-19. Additional rounds of docking, MD simulations and MM-PBSA calculations with remdesivir suggested that this compound can also work as a multi-target drug against SARS-CoV-2. Communicated by Ramaswamy H. Sarma

According to the literature SARS-CoV-2 belongs to the coronavirus family which has been around and causing pandemics for a while. However, this time the infection caused, called COVID-19, is a bit worse and capable of collapsing the lungs very quickly, killing the patient if no respiratory support is provided, being declared a worldwide pandemic in 11 March 2020 by the World Health Organization (WHO). In most countries the shortage of enough respirators, as well as, specialized personal to operate them, aroused a big fear of collapsing in the health systems. Therefore, social distancing and/or lockdown of the entire population was the most immediate response to this situation in order of smoothing the progress curve of the disease, and saving the Governments some time for preparing their health systems for facing properly this pandemic (Bedford et al., 2020; Cucinotta & Vanelli, 2020; Ortiz-Prado et al. 2020; Velavan & Meyer, 2020) .

However, improving the health infrastructure is not enough and the situation will be totally under control only after the naturally achieving of herd immunity by the population, the development of a vaccine, or emergence of a drug capable of avoiding the respiratory complications of COVID-19 (Fierabracci et al., 2020; Gates, 2020) . As the first two options can take months, maybe more than a year, the most immediate solution to beat COVID-19 is the search for new drugs based on the repurposing approach (Pushpakom et al., 2019) . On this sense we performed a comprehensive virtual screening (VS) study on a library comprising FDA approved drugs in a search for molecules capable of hitting the three main molecular targets of SARS-CoV-2 currently available in the protein data bank (PDB) (https://www.rcsb.org/) (Berman et al., 2000a) .

In order to select which proteins would be studied in terms of their interactions with small molecules, we analyzed the role that different proteins play in SARS-CoV-2 mechanisms of replication or cell entrance, and the availability and quality of their three-dimensional structures in the PDB (https://www.rcsb.org/) (Berman et al., 2000a) . All selected proteins are part of the viral machinery responsible for the replication and transcription of SARS-CoV-2 RNA genome or play important roles in the viral mechanisms of entrance inside human cells. Main protease (M pro ) is responsible for proteolysis of polyproteins 1a and 1ab, generating smaller proteins that are necessary for the virus replication Rathnayake et al., 2020) . Spyke protein is a homotrimer located in the external surface of the virus and contains the receptor-binding domain (RBD), which directly binds to the peptidase domain of the angiotensin-converting enzyme 2 (ACE2), and consequently promotes the entrance of the virus inside the cell (Cao et al., 2020; Yan et al., 2020) . Finally, papainlike protease (PL pro ) is responsible for processing the viral polyprotein, a process that is required for the release and maturation of several non-structural proteins, and thus is essential for the replication of SARS-CoV-2 and other coronaviruses (B aez-Santos et al., 2015; Gao et al., 2020; Klemm et al., 2020) .

Our search screened 1930 compounds contained in the FDA-approved drugs data set available at the Cheminformatic tools and Databases server (ChemoInfo) (https://chemoinfo. ipmc.cnrs.fr/) (Douguet, 2010 (Douguet, , 2018 against the targets cited above, and the 50 best molecules were selected for further docking studies using Molegro Virtual Docker (MVD V R ) 6.0 (Thomsen & Christensen, 2006) , in order to confirm their potential for drug repurposing against COVID-19. The 7 drugs that showed potential to bind the three protein targets selected in this study were submitted to molecular dynamics (MD) simulations and MM-PBSA calculations in order to corroborate the docking results. We also performed the same theoretical studies for remdesivir in complex with each target protein, since this drug has been reported in literature as presenting promising in silico and experimental results for the treatment of patients with COVID-19 (Grein et al., 2020; Hendaus, 2020; Naik et al., 2020; Wang et al., 2020) .

The molecular targets selected from the PDB (Berman et al., 2000b) and their respective binding pockets chosen to run the VS study, are listed on Table 1 . The search space of each molecular target in Table 1 were centered in the respective co-crystalized ligands.

The spherical search spaces for the VS studies on each target were determined with MVD V R (Thomsen & Christensen, 2006) as centered in the respective binding pockets defined in Table 1 . The coordinates and radius used for each search are listed in Table 2 .

As mentioned above we used the FDA-approved drugs data set (called here FDA-library) available at ChemoInfo server (https://chemoinfo.ipmc.cnrs.fr/) that provided 1930 FDAapproved drugs ready for docking (Douguet, 2010 (Douguet, , 2018 . Those molecules were evaluated regarding their interactions with each of the protein targets listed in Table 1 .

A Receptor-based Virtual Screening (RBVS) was carried out for each target protein separately using the ChemoInfo VS tool, which employs the docking algorithm Protein-Ligand ANT System (PLANTS) (Korb et al., 2009) . For each target protein, the respective .pdb file and coordinates of the center of the binding site, as defined in Table 2 , were given as inputs. After evaluation of the FDA-library, the ChemoInfo VS tool returned the results consisting on a list of the molecules ranked by the PLANTS (Korb et al., 2009 ) scoring function. The top results correspond to molecules that stablished the most favorable interactions with the target.

For each target, the best 50 results returned by the ChemoInfo VS tool were retrieved and these results were filtered in order to remove drugs that mainly act in the central and/or peripheral nervous system, renal system and/or cardiovascular system, because these effects would make it impossible to use those types of drugs in patients infected with COVID-19. After the filtering, we looked for drugs that showed potential to interact with more than one target protein and, therefore, triggering a synergistic effect, enhancing their potential to stop or slow down the SARS-CoV-2 action inside the human body. The drugs were then reranked in decreasing order of the number of protein targets that they interacted favorably with and, the ones showing favorable interactons with the three targets, were submitted to MVD V R and rounds of MD simulation and MM-PBSA calculations. A flowchart summarizing the VS procedure used in this study is illustrated in Figure 1 .

The drugs that showed favorable interactions with all three targets selected in this work were submitted to MVD V R , using the same protocol validated before (Botelho et al., 2020) , in order to re-evaluate the protein-ligand interactions using the Moldock Score, a docking algorithm more accurate than PLANTS (Korb et al., 2009) . For this task their 3 D structures were constructed in the Spartan 08 software (Hehre et al., 2006) , using the semi-empirical method Parametric Method 3 (PM3) (Stewart, 2004) for geometry optimization, and the Natural Population Analysis (NPA) (Reed et al., 1985) method for atomic charges calculations. For each protein-ligand complex, 10 runs were carried out in the search spaces defined in Table 2 , with 30 poses returned for each run. Those poses were analyzed regarding their MolDock Score and the residues they formed hydrogen bonds (H-bonds) with. After, the best pose of each drug in each target was selected to be submitted to MD simulations and MM-PBSA calculations. The complexes remdesivir/targets were also generated with MVD V R , following the same protocol described above and the best poses also selected for the further MD simulations and MM-PBSA calculations. The interaction energies of the ligands found inside the crystallographic structures listed in Table 1 were also determined using MVD V R for reference purposes.

The molecular dynamics (MD) simulations studies were performed applying the bonded and non-bonded parameters for the all-atom force field OPLS-AA (Kaminski et al., 2001) . The three-dimensional coordinates and topologies of the targets listed in Table 1 were generated by the pdb2gmx software, which is part of the GROMACS 5.1.4 package (Berendsen et al., 1995; Hess et al., 2008; Pronk et al., 2013; Van Der Spoel et al., 2005) , used in this work to perform all MD simulations. Since the OPLS-AA force field (Kaminski et al., 2001) has no parameters for the ligands studied, the AcPype software (Sousa Da Silva & Vranken, 2012), was used to generate topology and coordinates files of all ligands submitted to MD simulations. The Restrained Electrostatic Potential (RESP) method Cornell et al., 1993; Ribeiro et al., 2008; Wang et al., 2000) parameterized to reproduce the method Hartree-Fock and the basis set 6-31 G Ã (Jakalian et al., 2002) , was applied to calculate the atomic charges.

The target proteins complexed with the selected ligands were confined inside boxes under periodic boundary conditions (Mart ınez et al., 2007) where the effect of the solvent was reproduced by the water model TIP4P (Harrach & Drossel, 2014; Jorgensen et al., 1983; Mark & Nilsson, 2001) . For each target protein, the box format and size were selected seeking optimization of simulation time and absence of undesired artifacts. The boxes volumes and approximate numbers of water molecules of each system are listed in Table 3 . All systems were submitted to two minimization steps using the steepest descent algorithm with the convergence criterion of 100.00 kJ mol À1 nm À1 . The first step with position restrained (PR) and the second without PR. The equilibration of pressure and temperature was achieved through 100 ps of simulation using the canonical ensemble (NVT) (Bosko et al., 2005) , keeping the number of particles, volume and temperature constant. This was followed by 100 ps of simulation with an isothermal-isobaric ensemble (NPT) (Bosko et al., 2005) , keeping the number of particles, pressure, and temperature constant. After the equilibration step, all systems were submitted to a production step of 50 ns. All MD simulations were performed at 310 K and 1 bar, using 2 fs of integration time with the lists of pairs being updated at every 5 steps. The cut-off for Lennard Jones and Coulomb interactions were between 0 and 1.2 nm. The leap- frog algorithm was used in the production step with the Nose-Hoover thermostat (Evans & Holian, 1985) (s ¼ 0.5 ps) at 310 K and the Parrinello-Rahman barostat (Parrinello & Rahman, 1981 ) (s ¼ 2.0 ps) at 1 bar. All Arg and Lys residues were assigned with positive charges and all Glu and Asp residues were assigned with negative charges. The Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) software was used to visualize the simulation trajectories, and the Grace program (Turner, 2005) was used to build graphs of energy and root mean square deviation (RMSD). Additionally, the same procedure described above was used to run 100 ns of MD simulations with the three crystallographic complexes listed in Table 3 , in order to validate our MD simulation protocol.

MM-PBSA calculations of the binding free energies of each ligand after the MD simulations were carried out in order to support the former results. Through this method it is possible to predict the binding free energy of the ligands considering the vacuum potential energy, which includes both bonded and nonbonded interactions, as well as the free energy of solvation, which considers both polar and nonpolar terms (Kumari et al., 2014) . After MD simulations, the ligands were submitted to the g_mmpbsa tool (Kumari et al., 2014) 

The search for potential multitarget drugs returned the 7 drugs listed in Table 4 , which have the potential of interacting with all targets listed in Table 1 , suggesting that further repurposing studies of those drugs might lead to the discovery of an existing drug with potential to treat COVID-19. The structure of remdesivir is also shown in Table 4 . Some of the drugs presented in Table 4 have very similar structures that justify their similar results. Cobicistat and ritonavir, two antiviral drugs, share several common characteristics as well as nystatin and amphotericin B, which are also very similar to each other. Although having different structures, deferoxamine, itraconazole and elbasvir also showed potential to form stable complexes with all targets. Despite the structural differences, all drugs in Table 4 share the common feature of having several highly electronegative atoms linked to hydrogen atoms, i.e. spots where H-bonds can be formed, which probably contributed to their ability to form stable complexes with multiple proteins.

As mentioned above the 8 drugs presented in Table 4 were submitted to more docking calculations, with each molecular target chosen in this work, using MVD V R , so multiple poses per drug per protein could be analyzed regarding their stability.

The results of docking calculations using MVD V R shown in Figures 2-5 , Table S1 and Figures S1-S25 , confirmed the RBVS results, showing that the best poses of all 7 selected drugs are capable of binding to the targets, presenting large negative values of interacting energy (most below À100.00 kcal mol À1 ) comparable with the results for the reference ligands inside SC2M pro and SC2PL pro . In the case SC2Spyke the interaction energy of the reference ligand is not comparable to our ligands once the peptide LCB1 is a much larger ligand. This peptide is capable of establishing many more interactions with the RDB domain of SC2Spyke (see Figure S25 ), reflecting in the much more negative binding energy value shown in Figure 5 . Our ligands, however, presented binding energies with SC2Spyke similar to the other targets which is enough to qualify them as potential good binders to this protein. It was also observed that remdesivir showed similar energy values with all targets. This suggests that it can also qualify as a multi-target drug against SARS-CoV-2. Among the residues observed in interactions with the best poses of each ligand, it's possible to see in the 3D plots of Figures 2-4 and the 2D plots of Figures S1-S24 that most of the residues listed in Table 1 present some kind of interaction with the ligands and, also, that at least one residue of the respective binding site forms H-Bond with the ligand. This confirms that the selected drugs are capable of binding in the binding sites pointed in Table 1 and, therefore, have potential to inhibit the selected targets. Elbasvir was the best binder to all targets according to the docking results, with energy values always below À150 Kcal/mol, while deferoxamine presented the worse (highest) energy values with all targets except SC2MPro, where it was the second worse (Table S1 ). The other 5 ligands presented oscillating results among the targets, ranking well with some and not so well with others.

The selection of the best pose of each drug in each target resulted in 24 systems ligand/target that were submitted to additional rounds of MD simulations and MM-PBSA calculations.

In order to corroborate the docking results, 50 ns of MD simulations were performed with the 24 systems described above (one pose per drug of Table 4 in complex with each protein described in Table 1 ), plus the three crystallographic complexes listed in Table 3 , following the protocol described in the methodology section.

Plots of total and average energy of the crystallographic complexes ( Figure S26) show that the systems achieve stability since the beginning of the simulation and allowed reducing the simulated time to 50 ns in order to save computing resources.

Figures S27-S29 and the RMSD of each protein/ligand complex (Figure 6 ), confirm the stability of the systems during all simulated time. It's possible to see in those figures that the majority of the systems stabilized since the beginning of the MD simulations and that the RMSD of the ligands and proteins never oscillated over 0.7 nm. Deferoxamine and itraconazole were the ligands with worst behavior, showing larger RMSD oscillations with SC2Spyke (deferoxamine) and SC2M pro (itraconazole). Elbasvir (shown in blue in Figure 6 ) behaved well when in complex with all three proteins, indicating the favorable interactions between this ligand and the SARS-CoV-2 targets. Nystatin, shown in pink, presented stable RMSD values specially when complexed with SC2PL pro and SC2Spyke. As for SC2M pro , ritonavir (shown in purple) had a good behavior regarding RMSD values; pointing to the possibility of the administration of combined drugs for treating patients with COVID-19. Additionally, remdesivir, which has already shown promising clinical results for COVID-19 treatment (Grein et al., 2020; Hendaus, 2020; Naik et al., 2020; Wang et al., 2020) , seems to form more stable complexes with SC2PL pro and SC2Spyke than with SC2M pro , enlightening aspects of this drug action regarding SARS-CoV-2. The protein presenting worst behavior was SC2M pro when complexed with amphotericin B.

The H-bonds formed between each protein and each drug during the MD simulations were also analyzed in order to verify if such interactions, predicted by docking studies (Figures 2-4) , could also be observed during the MD simulations. The reference ligands were also analyzed in this aspect so the drugs could be compared to them in terms of Hbonds. Figures 7-9 show the H-bonds formed between each molecule and each protein. According to Figure 7 , the reference ligand to SC2M pro tends to form H-bonds mainly with Glu166 (pink) and Gln189 (green). Amphotericin B, deferoxamine and remdesivir seems to be more promising ligands than the reference in terms of H-bonds, since they interacted with a larger variety of residues during the MD simulations. Regarding interactions with SC2PL pro (Figure 8) , elbasvir is not shown because it did not form H-bonds during most of the MD simulation, although this drug presented good docking and RMSD results. As for the other ligands, none seem to be able to form numerous H-bonds with this protein; indicating that there may have non-polar interactions contributing for the complexes stabilities. Remdesivir seems to be the only drug capable of forming more H-bonds with SC2PL pro than the reference ligand. This results contributes to explain the already observed good clinical results of this drug. Considering SC2Spyke, Figure 9 shows that most drugs tend to form H-bonds mainly with Tyr449 (red) and Gln493 (blue), while the reference ligand formed H-bonds with multiple residues over the MD simulation time, mainly Lys417 (cyan), Asn487 (dark green), Gly502 (yellow), and Tyr449 (red). This is not unexpected since the reference ligand is a peptide and has nearly 60 amino acids (Cao et al., 2020) , so its size favors the formation of multiple interactions with SC2Spyke residues. Among the studied drugs, deferoxamine, nystatin and remdesivir are the ones that formed H-bonds with SC2Spyke residues during most of the MD simulation time, indicating that they tend to stay bound to this protein. 

The results of the MM-PBSA calculations are shown in Figure  10 . The negative binding energies predicted for all ligands (most below À50 Kcal mol À1 ) with all targets corroborated the docking and MD simulations results, confirming their stability and capacity of binding to these targets. Reference ligands of SC2M pro and of SC2PL pro had similar results to the other drugs, confirming their potential to inhibit those proteins. SC2Spyke reference ligand presented a significantly lower binding energy than all the other ligands, probably due to its high similarity with the SC2Spyke natural substrate (Cao et al., 2020) and consequent extremely low IC 50 value (Table 3) . However, other drugs such as elbasvir, itraconazole, nystatin and ritonavir also had good MM-PBSA results for SC2Spyke. Comparing the drug performances for the three proteins, elbasvir and ritonavir are the drugs presenting best general results.

We performed a comprehensive theoretical study applying the methods of RBVS, docking, MD simulations and MM-PBSA calculations, that allowed the selection of 7 potential hits amongst 1930 drugs approved by the FDA (https:// chemoinfo.ipmc.cnrs.fr/MOLDB/index.php#ref1). Our 7 hits, together with remdesivir, showed capacity of effectively binding to the three main SARS-CoV2 targets SC2M pro , SC2PL pro and SC2Spyke, qualifying as potential preapproved multi-target drugs. These results suggest that these drugs, worth being immediately submitted to in vitro evaluation against SARS-CoV2. Also, the fact of being multi-target hits suggests that these compounds could have a synergistic effect if administered simultaneously against COVID-19. The most promising drugs for this, according to our results, are remdesivir, elbasvir and ritonavir.

No potential conflict of interest was reported by the authors. 

Productivity Researcher Program PPP); and FAPES, grant number 03/ 2020-2020-WMT5F. This work was also supported by the University of Hradec Kralove (VT2019-2021). We also thanks the Federal University of Lavras (UFLA) for software facilities

Kamil Kuca

The SARS-coronavirus papain-like protease: Structure, function and inhibition by designed antiviral compounds

A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model

COVID-19: towards controlling of a pandemic

GROMACS: A message-passing parallel molecular dynamics implementation

The protein data bank

The protein data bank

Molecular simulation of dendrimers and their mixtures under shear: Comparison of isothermal-isobaric (NpT) and isothermal-isochoric (NVT) ensemble systems

Ligand-based virtual screening, molecular docking, molecular dynamics, and MM-PBSA calculations towards the identification of potential novel ricin inhibitors

De novo design of picomolar SARS-CoV-2 miniprotein inhibitors

Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation

WHO declares COVID-19 a pandemic

Molecular modeling studies on the interactions of 7-methoxytacrine-4-pyridinealdoxime, 4-PA, 2-PAM, and obidoxime with VX-inhibited human acetylcholinesterase: A near attack conformation approach

Molecular modeling studies on the interactions of aflatoxin B1 and its metabolites with the peripheral anionic site of human acetylcholinesterase

E-LEA3D: A computational-aided drug design web server

Data sets representative of the structures and experimental properties of FDA-approved drugs

The nose-hoover thermostat

COVID-19: A review on diagnosis, treatment, and prophylaxis

Crystal Structure of SARS-CoV-2 papain-like protease

Responding to Covid-19 -A once-in-a-century pandemic?

Compassionate use of remdesivir for patients with severe Covid-19

Structure and dynamics of TIP3P, TIP4P, and TIP5P water near smooth and atomistic walls of different hydroaffinity

Remdesivir in the treatment of Coronavirus Disease 2019 (COVID-19): A simplified summary

GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation

VDM: Visual molecular dynamics

Fast, efficient generation of high-quality atomic charges. AM1-BCC Model: II. Parameterization and validation

Structure of Mpro from COVID-19 virus and discovery of its inhibitors

Comparison of simple potential functions for simulating liquid water

Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides

Mechanism and inhibition of the papain-like protease, PLpro, of SARS-CoV-2

Empirical scoring functions for advanced protein-ligand docking with PLANTS

g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K

Fundamentos de Simulac¸ão Por Dinâmica Molecular

Remdesivir (GS-5734) as a therapeutic option of 2019-NCOV main protease -in silico approach

Clinical, molecular, and epidemiological characterization of the SARS-CoV-2 virus and the Coronavirus disease 2019 (COVID-19), a comprehensive literature review

Polymorphic transitions in single crystals: A new molecular dynamics method

GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit

Drug repurposing: Progress, challenges and recommendations

3C-like protease inhibitors block coronavirus replication in vitro and improve survival in MERS-CoV-infected mice

Natural population analysis

MKTOP: A program for automatic construction of molecular topologies

ACPYPE -Antechamber python parser interface

Optimization of parameters for semiempirical methods IV: Extension of MNDO, AM1, and PM3 to more main group elements

MolDock: A new technique for high-accuracy molecular docking

XMGRACE, version 5.1.25

GROMACS: Fast, flexible, and free

The COVID-19 epidemic

How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules

Remdesivir in adults with severe COVID-19: A randomised, double-blind, placebo-controlled, multicentre trial

Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2