key: cord-1008860-s01yopm5
authors: Nunes, Vinicius S.; Paschoal, Diego F. S.; Costa, Luiz Antônio S.; Santos, Hélio F. Dos
title: Antivirals virtual screening to SARS-CoV-2 non-structural proteins
date: 2021-05-05
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2021.1921033
sha: 8cb25aad14df16b57094006901511bf5ea86fc18
doc_id: 1008860
cord_uid: s01yopm5

In March 2020, the World Health Organization (WHO) declared coronavirus disease-19 (COVID-19), caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a pandemic. Since then, the search for a vaccine or drug for COVID-19 treatment has started worldwide. In this regard, a fast approach is the repurposing of drugs, primarily antiviral drugs. Herein, we performed a virtual screening using 22 antiviral drugs retrieved from the DrugBank repository, azithromycin (antibiotic), ivermectin (antinematode), and seven non-structural proteins (Nsps) of SARS-CoV-2, which are considered important targets for drugs, via molecular docking and molecular dynamics simulations. Drug–receptor binding energy was employed as the main descriptor. Based on the results, paritaprevir was predicted as a promising multi-target drug that favorably bound to all tested Nsps, mainly adipose differentiation-related protein (ADRP) (-36.2 kcal mol(−1)) and coronavirus main proteinase (Mpro) (-32.2 kcal mol(−1)). Moreover, the results suggest that simeprevir is a strong inhibitor of Mpro (-37.2 kcal mol(−1)), which is an interesting finding because Mpro plays an important role in viral replication. In addition to drug–receptor affinity, hot spot residues were characterized to facilitate the design of new drug derivatives with improved biological responses.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a new beta coronavirus responsible for the worldwide pandemic of coronavirus disease-19 (COVID-19) since December 2019. Although it is similar to SARS-CoV ($80% sequence identity), SARS-CoV-2 has some subtle mutations that cause the highly contagious and severe COVID-19 (Walls et al., 2020) . Most infected people have mild clinical symptoms; however, they still need medical care, which can cause the healthcare system to collapse. Based on previous pandemic experiences and the recommendation of the World Health Organization (WHO), initially, some nonpharmacological interventions were established, such as social distancing, closing schools, and wearing face masks (Cobey, 2020; Leung et al., 2020; Liu et al., 2020; Verity et al., 2020) . Moreover, the scientific community has been dedicated to testing existing drugs (drug repositioning) for the treatment of COVID-19 and elucidating most of the molecular components of the virus structure. For drug repositioning, three main drugs have been considered: chloroquine, an antimalarial agent ; remdesivir, currently in clinical development to treat Ebola infection (Wang, 2020) ; and dexamethasone, a corticosteroid that reduced mortality by approximately one-third for patients on ventilators (Gordon et al., 2020) .

SARS-CoV-2 genome has a length of 30 kb and is very similar to the SARS-CoV genome. They have an open reading frame 1ab (Orf1ab) encoding 16 predicted non-structural proteins (Nsps) and four typical coronavirus structural proteins. The 5 0 Orf1ab gene encodes polyproteins that are auto-proteolytically processed into 16 Nsps, forming a replicase-transcriptase complex (RTC) (Gordon et al., 2020; Wu et al., 2020) . The RTC consists of multiple enzymes: the complex papain-like protease (PLpro), adenosine diphosphate ribose-1"-monophosphatase (ADRP) (Nsp3), main protease (Mpro) (Nsp5), Nsp7-Nsp8 primase complex, ribonucleic acid (RNA)-replicase (Nsp9), primary RNA-dependent RNA polymerase (RdRp) (Nsp12), helicase-triphosphatase (Nsp13), exoribonuclease (Nsp14), endoribonuclease (Nsp15), and N7/ 2'O-methyltransferases (Nsp10/Nsp16). The 3 0 end of the viral genome has four structural proteins (namely, spike (S), envelope (E), membrane (M), and nucleocapsid (N)) and nine putative accessory factors Pant et al., 2020) . High-resolution crystal structures of the many Nsps of SARS-CoV-2 have recently been reported in the Protein Data Bank (PDB).

Currently, several molecular docking and molecular dynamics (MD) simulation studies on different Nsps and drugs are being reported every week worldwide. Most of these studies involve repositioning of drugs, mainly antiviral drugs (Alamri et al., 2020; Elfiky, 2020; Gordy et al., 2020; Hakmi, 2020; Pant et al., 2020; Peele et al., 2020; Shah et al., 2020; Tazikeh-Lemeski et al., 2020; Wu et al., 2020) . These studies have demonstrated the inhibition potential of some antiviral drugs against SARS-CoV-2 Nsps. Recently, a study of infected patients at a Hong Kong hospital showed that a cocktail of three antiviral drugs (namely, interferon beta-1b, lopinavir-ritonavir, and ribavirin) was effective in treating patients infected with SARS-CoV-2 (Hung et al., 2020) . Therefore, antiviral drugs are becoming increasingly important in the repositioning of drugs for the treatment of COVID-19.

Herein, we implemented a structure-based virtual screening approach to repurpose available antiviral drugs as potential treatments for COVID-19. Drug repositioning is a proven, cost-effective, and time-saving solution for finding new indications for already established drugs (Hakmi, 2020; Peele et al., 2020; Shah et al., 2020; Tazikeh-Lemeski et al., 2020) .

Our strategy focuses on the identification of small molecules with potential activity against SARS-CoV-2 Nsps from a subset of protease inhibitors intended for the treatment of other viral infections, including hepatitis C virus (HCV) or human immunodeficiency virus (HIV).

Crystal structures of PLpro (PDB: 6W9C), ADRP (PDB: 6W02), Mpro (PDB: 6W63), RNA-replicase Nsp9 (PDB: 6W9Q), RdRp (Nsp12 -PDB: 6M71), endoribonuclease Nsp15 (PDB: 6WLC), and 2'O-methyltransferase Nsp16 (PDB: 6W4H) were selected and downloaded from the PDB. As some fragments of residues were missing in the 6W9C and 6M71 structures, it was necessary to generate homology models for these structures. The models were constructed using the Modeler software v9.23 (Eswar et al., 2007) . Subsequently, all seven structures were prepared and converted to Protein Data Bank, Partial Charge (Q), & Atom Type (T) (PDBQT) format files using AutoDockTools v1.5.7 (Morris et al., 2009) .

Nsp3 has 1945 residues that are divided into two domains: PLpro and ADRP (Lei et al., 2018; Yoshimoto, 2020) . PLpro is responsible for cleavages located at the N-terminus of the replicase polyprotein. In addition, PLpro has a deubiquitinating/deIS gylating activity and processes both Lys48and Lys63-linked polyubiquitin chains from cellular substrates (Angeletti et al., 2020; Yoshimoto, 2020) . It also participates in the assembly of virally induced cytoplasmic double-membrane vesicles necessary for viral replication. PLpro antagonizes innate immune induction of type I interferon by blocking the phosphorylation, dimerization, and subsequent nuclear translocation of host interferon regulatory factor-3 and prevents host nuclear factor-kappa-B signaling (Lei et al., 2018; Yoshimoto, 2020) . Nsp3 also exhibits ADRP activity, which has been proposed to play a regulatory role in the replication process (Xu et al., 2009; Yoshimoto, 2020) .

Nsp5 or Mpro has 306 residues and cleaves at 11 distinct sites to yield mature proteins of SARS-CoV-2 (Yoshimoto, 2020) . Mpro recognizes the sequence of Leu-Gln cleavage sites. Inhibiting the activity of this enzyme would block viral replication .

Nsp9, with 113 residues, is one of the smallest SARS-CoV-2 Nsps (Yoshimoto, 2020) . It participates in viral replication by acting as a soluble RNA-binding protein (Krichel et al., 2020; Littler et al., 2020) .

Nsp12 or RdRp has 932 residues and is the second largest Nsp (Yoshimoto, 2020) . It is responsible for the replication and transcription of the viral RNA genome (Peng et al., 2020; Yin et al., 2020) . Although Nsp12 is capable of conducting polymerase reactions, it demonstrates extremely low efficiency for these reactions. In contrast, the presence of Nsp7 and Nsp8 remarkably stimulates the polymerase activity of Nsp12 (Peng et al., 2020; Subissi et al., 2014) . Thus, the Nsp12-Nsp7-Nsp8 subcomplex is defined as the minimal core component for mediating coronavirus RNA synthesis (Peng et al., 2020) .

Nsp15 has 346 residues and shows endoribonuclease activity (Yoshimoto, 2020) . It is an Mn þ2 -dependent and uridylate-specific enzyme, which leaves 2 0 -3 0 -cyclic phosphates 5 0 to the cleaved bond Yoshimoto, 2020) .

Nsp16 has 298 residues and is an RNA-methyltransferase (MTase) that mediates mRNA cap 2 0 -O-ribose methylation to the 5 0 -cap structure of viral mRNAs. An N7-methyl guanosine cap is a prerequisite for the binding of Nsp16 (Tazikeh-Lemeski et al., 2020; Yoshimoto, 2020) and, thus, plays an essential role in viral mRNA cap methylation, which is essential for evading the immune system (Tazikeh-Lemeski et al., 2020) .

A set of 22 antiviral drugs, azithromycin (antibiotic), and ivermectin (antinematode) were selected from the DrugBank database indicated for COVID-19 ( Figure S1 ). Then, 12 structures were downloaded from the Cambridge Crystallographic Data Centre database, and 12 structures were created using the information obtained from the DrugBank database. The structures were optimized and characterized as a minimum point on the potential energy surface in the gas phase at the semi-empirical PM3 level employing the GAUSSIAN 09 Rev. D.01 program (Frisch et al., 2009 ).

Molecular docking simulations of all seven protein receptors against the approved set of 24 drugs were performed using the latest version of Autodock Vina v1.1.2 software (Trott & Olson, 2009 ). The number of poses was fixed at 20, whereas the remaining parameters were kept at their default values (exhaustiveness ¼ 8). Only polar H atoms were considered for proteins with all residues in their standard protonation states at neutral pH. The grid sizes and center coordinates of the seven boxes are presented in Table S1 . All seven grid boxes were set at the active site of the proteins. The results of the virtual screening experiment were ranked according to the binding affinity (BA) of the best scoring conformation ( Table S2 ). The top 10 ranked candidates were selected for further analysis (Table S3 ).

To improve the accuracy of docking simulations, the complexes of the best three candidates, namely, paritaprevir (PAR), simeprevir (SIM), and glecaprevir (GLE), obtained via the docking simulations (bold entries in Table S3 ) and the seven targets were subjected to MD simulations (a total of 21 structures were simulated at this stage). The proteins were prepared using the pdb4amber utility in AMBER package (Case et al., 2016) considering that the titratable residues were in their protonation states at neutral pH, that is, deprotonated asparagine (Asn) (AS4 -pKa ¼ 4.0) and glutamic acid (Glu) (GL4 -pKa ¼ 4.4) and protonated cysteine (Cys) (CYS -pKa ¼ 8.5), tyrosine (Tyr) (TYR -pKa ¼ 9.6), histidine (His) (HIP -pKa ¼ 6.6), and lysine (Lys) (LYS -pKa ¼ 10.4).

The final PDB files contained all the residues sequentially numbered from 1 to N. Table S4 demonstrates the procedure to convert the abovementioned numbered residues into original residues. For ligands, we used a set of parameters from the generalized AMBER force field 2 (GAFF2) as implemented in the AMBER package and semiempirical AM1-BCC charges, which were parameterized to reproduce HF/6-31G Ã restrained electrostatic potential (RESP) charges (Case et al., 2016) . We used the force field ff14SB (Maier et al., 2015) for the targets and the transferable intermolecular potential model (TIP3P) (Jorgensen et al., 1983) for solvents. The system containing the ligand and target in distinct docking modes was neutralized and immersed in a truncated octahedral water box with the minimum distance between the solute and the periodic box set to 8.0 Å. The MD simulations involved the following steps: (i) optimization of the solvent while keeping the solute fixed, (ii) optimization of the entire system, (iii) NVT heating up to 310.0 K for 250 ps (s ¼ 2 fs), (iv) NPT density equilibration for 1 ns (s ¼ 2 fs), and (v) three independent NPT runs for 20 ns each (s ¼ 1 fs). Temperature and pressure were controlled using a Berendsen thermostat (Berendsen et al., 1984) and a Langevin barostat (Uberuaga et al., 2004) , respectively. The non-bonded cutoff was 12 Ð, and long-range electrostatic interactions were evaluated using the particle mesh Ewald approach (Ryckaert et al., 1977) . All bonds involving H atoms were constrained via the SHAKE algorithm (Ryckaert et al., 1977) . Production trajectories were analyzed and employed to predict the ligand-receptor binding energy using the generalized Born/surface area model (GB/SA) (Miller et al., 2012) with igb ¼ 2 in the sander program. The binding energy for each run was averaged over 200 frames evenly spaced in time.

Finally, the previously calculated GB/SA binding energies were used to narrow the screening and to select two best ligands and three best targets. Subsequently, the six complexes were subjected to a longer MD simulation using the abovementioned protocol with the production phase extended to 200 ns with s ¼ 2 fs. The last 100 ns were considered for analyzing the GB(PB)/SA binding energy (PB stands by the Poisson À Boltzmann approach), root-mean-square deviation (RMSD), frequency of contacts, and ligand-receptor H bonds.

In summary, the virtual screening was initiated with 24 ligands and seven targets, and the number of ligands and targets was progressively reduced to three and seven using the docking results, respectively. Then, based on the MD simulation results, the best two ligands and three targets were chosen, which were further examined via long MD simulations. The analyses of trajectories were performed using CPPTRAJ software (Roe & Cheatham, 2013) , and images created by the Visual Molecular Dynamics (Humphrey et al., 1996) and Discovery Studio Visualizer V.20.1.0.19295 programs.

A previously reported study on the SARS-CoV-2 genome (Lu et al., 2020) allows drug designers to identify the main viral proteins as potential targets and search for possible drugs to prevent the spread of the virus. The main targets are Nsps because they are directly involved in the replication and maturation of SARS-CoV-2 inside the host cell. Table S3 presents the results of the docking simulations for the 10 best candidates among the 24 drugs evaluated herein (see Table S2 for the entire set of ligands).

The next step of virtual screening was to select antiviral drugs with BA lower than À9.0 kcal mol À1 , and herein, PAR, SIM, and GLE were chosen (entries in bold in Table S3 ). PAR showed the best BA for all the seven selected targets. PAR is a part of combination therapy for the treatment of chronic hepatitis C, an infectious liver disease caused by infection with hepatitis C virus (HCV) (Badri et al., 2016) . Its action prevents viral replication by inhibiting the NS3/4A serine protease of HCV; PAR was first available in the market in a fixeddose combination with ombitasvir, dasabuvir, and ritonavir as an FDA-approved medication (Badri et al., 2016) . Its common side effects include headache, fatigue, nausea, pruritus, and insomnia, reported in >10% of patients (Kumar et al., 2019) . Other studies have reported a lack of significant side effects and a short duration of therapy, which are considerable advantages of PAR over other antiviral drugs (Alamri et al., 2020; Yaras et al., 2019) . A virtual screening study conducted using Autodock Vina software demonstrated that PAR was the best antiviral drug for inhibiting Mpro (Khan et al., 2020) , and the estimated BA was À9.8 kcal mol À1 .

Other antiviral drug selected from docking simulations was SIM. It showed the second-best BA to four of the seven selected targets, that is, PLpro, Nsp12, Nsp15, and Nsp16 (Table S3) . SIM is a NS3/4A serine protease inhibitor indicated in patients with HCV genotype 1 and has been classified as a second-generation protease inhibitor (Alamri et al., 2020; Mahdian et al., 2020). As the viral NS3/4A protease complex is essential for cleaving the HCV-encoded polyprotein into individual viral proteins facilitating replication, SIM blocks the viral replication process (Mahdian et al., 2020) .

The third antiviral drug selected in this study was GLE because it exhibited a suitable BA (-9.2 kcal mol À1 ) to ADRP (Table S3) . GLE is a direct-acting antiviral agent against HCV and an NS3/4A protease inhibitor that targets viral RNA replication. It is a useful drug for patients who experience therapeutic failure from other NS3/4A protease inhibitors (Lawitz et al., 2016) . Recently, via well-defined computational methods, the combination of GLE and maraviroc was identified as a potential inhibitor of SARS-CoV-2 Mpro (Shamsi et al., 2020) . Both drugs bind to the substrate-binding pocket of SARS-CoV-2 Mpro and form a significant number of noncovalent interactions with conserved residues (Shamsi et al., 2020) .

The 21 complexes (three ligands and seven receptors) selected via the docking simulations were subjected to MD simulations (60 ns) in an aqueous solution under thermodynamic conditions (T ¼ 310 K) to relax the docking poses and narrow the antiviral screening. Herein, three distinct trajectories of 20 ns each were analyzed for each drug-receptor complex, and a total of 6,000 frames were monitored over 60 ns. The average GB/SA binding energies for all complexes are reported in Table 1 and plotted against the Vina docking scores in Figure 1 .

Correlation between the data shown in Figure 1 suggests that the binding mode for most ligands is not substantially affected by dynamics in a short simulation time; however, some deviations were observed. Considering the most stable complexes, PAR was predicted as the strongest inhibitor for five of the seven receptors among the tested antiviral drugs (Table 1 ). In contrast, for Mpro and Nsp15, SIM and GLE were predicted to be the best ligands, respectively. With respect to receptors, PAR and SIM demonstrated the highest BAs for ADRP, Mpro, and Nsp12 and were the best candidates as inhibitors for these receptors, with the GB/SA binding energy < À40 kcal mol À1 . Therefore, the two best inhibitors PAR and SIM and the three most sensitive targets ADRP, Mpro, and Nsp12 were selected for further analysis and discussion based on long (200 ns) MD simulations.

The six drug-receptor complexes selected from the sequential molecular docking-MD simulations (bold entries in Table 1 ) were subjected to 200 ns of simulation. Structural evolution over time is represented by the RMSD in Figure 2 . The images presented in the figure were obtained at every 20 ps during a 100 ns simulation, with the first frame considered as a reference for each trajectory. Note that all structures stabilized after 20 ns (1000 frames) and the average RMSD was lower than 2.8 Å. Most of the structural changes originated from ligands. For instance, for PAR@Nsp12, which presented the largest average RMSD (2.8 Å), the RMSD values for the protein backbone (Ca only) and PAR were 2.2 and 2.9 Å, respectively. The complexes of drugs with Nsp12 were more flexible, followed by those with Mpro and ADRP (see the caption of Figure 2 for the corresponding average RMSDs). Subsequently, we investigated the effects of the high flexibility of the drug-Nsp12 complexes on the overall stability of these complexes. Average GB/SA binding energies calculated over the last 100 ns are presented in Table 2 , and the main individual contributions are explicitly included.

By comparing the values provided in Tables 1 and 2, we noticed that except for PAR@Mpro and SIM@ADRP, which were slightly more stable during the long simulation, the remaining four complexes were less stable at long simulation times. Moreover, the stability order of the complexes with ADRP and Mpro was the same as exhibited in Tables 1 and 2 and was different from that of the complexes with Nsp12, which were predicted to be over-stabilized in a short simulation time. For easy comparison, the data shown in Tables 1 and 2 are graphically presented in Figure S2 .

Considering the importance of Mpro in viral replication, high stability of SIM@Mpro is crucial. This complex was the most stable among all the complexes evaluated herein, and its binding energy in aqueous solution was À37.2 kcal mol À1 . Most of the binding energy (66%) is contributed by non-electrostatic interactions (van der Waals (v.d.W.), Table 2 ) because of the large size of the ligand and thus significant number of dispersive contacts. Nevertheless, electrostatic contribution plays a primary role in the relative stabilities of SIM@Mpro and PAR@Mpro, accounting for 34% of the total binding energy. The electrostatic contribution to the binding energy of SIM@Mpro was approximately twice that calculated for PAR@Mpro (Table 2 ). Figure 3 shows the GB/SA binding energy per frame for the four most stable complexes. The values were calculated for 250 frames sampled at 400 ps intervals along the trajectory. The binding energy of PAR@Mpro slightly decreased with time mainly because of the increase in the electrostatic contribution (Figure 3a and b). In contrast, SIM@Mpro stabilized over time owing to the decrease in the non-electrostatic contribution, mainly at the end of the simulation. Overall, Figure 3a and b clearly À53.5 ± 6.4(±0.4) À13.9 ± 9.4(±0.6) À67.5 ± 13.4(±0.8) À32.2 ± 5.7(±0.40) SIM@Mpro À51.1 ± 5.0(±0.32) À25.8 ± 6.4(±0.40) À76.9 ± 8.1(±0.50) À37.2 ± 4.6(±0.29) PAR@Nsp12

À39.6 ± 5.9(±0.37) À25.4 ± 10.6(±0.70) À64.9 ± 13.8(±0.90) À16.9 ± 4.4(±0.28) SIM@Nsp12

À30.2 ± 6.4(±0.40) À8.99 ± 6.4(±0.40) À39.2 ± 9.8(±0.62) À18.5 ± 5.3(±0.33)

The values are in kcal mol À1 . DE v.d.W and DE ele are the van der Waals and electrostatic contributions for the binding energy, respectively. DG gas is the binding free energy in gas phase (¼ DE v.d.W þ DE ele ) and DG aq is the binding free energy in aqueous solution (¼ DG gas þ DG solv ). The binding energy was calculated over 250 frames. The values in parenthesis are the standard mean error (SME) and reflect how precise the mean value is as an estimate of the true mean.

demonstrate that although the non-electrostatic contacts drive the binding of both antiviral drugs to the receptors, electrostatic contacts are essential for the binding of SIM to Mpro. To gain insights into the molecular bases of the binding mode and complex stability, the ligand-receptor contacts were mapped over the MD trajectory. The native contacts were labeled as any two atoms within a cut-off distance of 5 Å in a reference structure, which was defined as the first frame of the last 100 ns of the production trajectory. The results for SIM@Mpro and PAR@Mpro are shown in Figure 4a and c, respectively, where the total fraction of frames in which the specified contact is present (frequency of contacts) is plotted against the number of residues. Figure 4a and c might be used to establish the binding sites. Further, in this specific case, the binding sites were represented by the residues with a frequency of contacts > 20%. Note that the frequency of contacts is based on geometric criteria, and, hence, the contacts can be either attractive or repulsive. Thus, the frequency of contacts provides an idea of the "strength" of the ligand-receptor interaction.

In the most stable complex SIM@Mpro (Figure 4a) , the important ligand-receptor contacts were in Asn142-Cys145, Glu166, and Gln189 enzymatic moieties. In Asn142-Cys145, some contacts are assigned as H bonds (Figure 4b ). H bond was defined based on the structural criteria of a donor (D-H)-acceptor (A) distance of 3.0 Å and D-H-A angle cutoff of 135 . Figure 4b shows the frequency of the H bonds formed during the trajectory per residue. The three main anchor points, which may be considered the minimum requirement for drug activity, are Gly143 (59%), Cys145 (22%), and Asn142 (12%), with SIM acting as a H-bond acceptor. Moreover, the H bond with Gly143 has the longest lifetime (21); that is, it appears in 21 consecutive frames (420 ps). A representative structure of SIM@Mpro is shown in Figure 5a and b, where the residues participating in H bond formation are indicated.

Asn142-Gly143-Ser144-Cys145 forms a network of H bonds around the cyclopropylsulfonyl group of SIM, stabilizing the SIM-receptor complex. On the other side of the binding site, the hydrophobic thiazole-quinoline group is in frequent contact with Gln189, strengthening the v.d.W. interaction. The thiazole-quinoline group is exposed to the solvent (solvent-accessible surface area (SASA) ¼ 416.9 Å 2 ) and has considerable mobility along the MD trajectory. For PAR@Mpro, we observed large movements in the methylpyrazine group. However, the phenanthridine group remained at the same binding site along the entire trajectory ( Figure  5c) , which was the main anchor point for PAR@Mpro docking. The H bond with Gln189 appears in 24% of the frames (Figure 4d ) with a short lifetime (<14 consecutive frames, 280 ps). Thr24 and Ser46 also form low-frequency H bonds with the methylpyrazine group of PAR; nevertheless, these bonds have a short lifetime due to the high flexibility of the methylpyrazine group. Consequently, 79% of the favorable binding energy of PAR@Mpro originates from v.d.W. according to Table 2 .

Note that SIM and PAR bind to Mpro in different ways ( Figure 5 ). Although both ligands have a cyclopropylsulfonyl group, which is strongly H bonded to Asn142-Gly143-Ser144- Cys145 in SIM@Mpro, the docking mode in PAR@Mpro is driven by the intercalation of the phenanthridine group into the protein cleft between Leu141-Asn142 and Glu166 residues. Accordingly, two main molecular characteristics can explain this difference: the methylpyrazine group in PAR favorably interacts with Thr24 and Ser46, positioning the hydrophobic phenanthridine group toward Asn142-Gly143-Ser144-Cys145, and the large volume of the thiazole-quinoline group in SIM prevents intercalation of this group into the protein gap formed by the Leu141-Asn142 and Glu166 residues. Moreover, the energy gain in PAR@Mpro due to v.d.W. interactions is the same as that predicted for SIM@Mpro ($-50 kcal mol À1 , Table  2 ). This implies that the binding mode of SIM must be predominantly driven by the H bonds between the cyclopropylsulfonyl group and the Asn142-Gly143-Cys145 triad. To obtain a comprehensive understanding, a per-residue energy decomposition analysis was carried out. Figure 6 shows the GB/SA energy contributions for the 10 main residues. The hot spot residues (with a binding energy contribution of < À1 kcal mol À1 ) are shown in the red bar.

Results for Mpro are shown in Figure 6a and b. For SIM@Mpro (Figure 6b) , the hot spot residues Asn142-Gly143-Ser144-Cys145 and Met49 contribute 28% of the total binding energy. Met49 attracted our attention as it had a significant contribution to the binding energy owing to the short contact of its sidechain with the macrocycle ring in the ligand. The energy contribution of Met49 completely arises from the v.d.W. term, whereas the energy contribution of the Asn142-Gly143-Ser144-Cys145 residues, which form H bonds with the ligand, originates from the Coulomb term. For PAR@Mpro (Figure 6a ), Asn142 and Gln189 are the main anchor points; Asn142 is responsible for positioning the phenanthridine group (Figure 5c) , and, therefore, its contribution to the binding energy stems from the non-electrostatic term. Meanwhile, the Gln189 residue is a part of the H bond formed with the ligand (Figure 4d) , and, thus, its contribution to the binding energy originates from the electrostatic term. Similar to the case of SIM@Mpro, Met49 acts as an anchor point for PAR@Mpro binding.

Mpro is certainly the main drug target in SARS-CoV-2 because its inhibition blocks viral replication. In a recent study, Kumar et al. (Kumar et al., 2020) evaluated 13 antiviral drugs against Mpro using in silico approaches. They found that indinavir was the best inhibitor, with a docking score of À8.8 kcal mol À1 , among the analyzed antiviral drugs. Its mode of interaction involves the formation of H bonds with Glu166 and Gln189, which are also found in frequent contacts with SIM and PAR (Figure 4 ). In addition, in the same study, a broad set of hydroxyethylamine (HEA) analogs was examined; the main compound (labeled as 16) formed H bonds with Gly143, Ser144, Cys145, and His164, and the H bond with Gly143 was maintained throughout the MD trajectory. According to the present study, the H bonds of SIM with Gly143 and Cys145 were retained for more than 20% of the simulation time (Figure 4b) , and the H bond of SIM with Gly143 was observed in 60% of the 5,000 frames analyzed along 100 ns of simulation. Besides, Figure 6b shows these two residues are the main hot spots in the binding mode. Figure 6 . Per-residue energy decomposition analysis. The so-called hot spot residues (contributing to the binding energy by < -1 kcal mol À1 ) are shown in red bar.

Interestingly, SIM interacts with Mpro via the main anchor points found for indinavir and HEA 16 (Shamsi et al., 2020) ; this can help understand the most favorable binding energy. Previously, hierarchical virtual screening was used to test more than 2,000 approved drugs from DrugBank against Mpro, among which 39 drugs were further subjected to MD simulations (Wang, 2020) . Results suggested that Asn142, Glu166, and Gln189 were hot spot residues. Herein, in addition to Leu141, Asn142, and Glu166, Gln189 was assigned as an important anchor point for PAR@Mpro binding (Figure 4d and 6a). Huynh et al. (Huynh et al., 2020) used in silico methods to evaluate 19 market drugs against Mpro and chose entecavir, used to treat hepatitis B virus infection, for MD simulation. Their findings showed that the main anchor points were His41, Glu166, and Gln189, as found for PAR in the present study (Figures 5c and d & 6a) . Recently, Gao et al. reported a very interesting study based on the repositioning of structure-based drugs to more than 8,000 molecules using Mpro as a target. Proflavin was the best representative of the FDA-approved set of drugs and binds Mpro via Glu166 as an anchor point. The other two promising candidates are chloroxine and demexiptyline, which also bind to Mpro via Ser144 and Cys145. These results indicate that SIM and PAR possess the main structural requirements for strongly binding to Mpro, as found in the present study. Moreover, the results acquired for SIM@Mpro suggest that the orientation of the polar group acting as a H-bond acceptor close to the Gly143-Ser144-Cys145 triad might be a relevant molecular feature to the design of novel Mpro inhibitors (Figure 5a and b) . Although PAR also has a polar cyclopropylsulfonyl group in its structure as in SIM, its binding mode is driven by the hydrophobic methylpyrazine and phenanthridine groups (Figure 5c and d) .

Stability of the complexes with ADRP was slightly lower than that of the complexes with Mpro, and PAR@ADRP was the second most stable complex ( Table 2 ). The GB/SA binding energies for SIM@ADRP and PAR@ADRP were À27.6 and À36.2 kcal mol À1 , respectively. Higher stability of PAR@ADRP is mainly because of the dispersion contribution, which account for 85% of the favorable binding energy (Figure 3c ). This favorable hydrophobic contact between PAR and ADRP is due to the interaction of the phenanthridine group of PAR with Ala127, Phe130, and Pro154 ( Figure 7c ). Electrostatic contributions for PAR@ADRP are only À10.6 kcal mol À1 because of the formation of H bonds between Phe154 and the cyclopropylsulfonyl group (24%) and Gly46 (12%) and between Leu124 and the methylpyrazine-2-carbonyl group (5%) (Figure 7d ). Although these H bonds have measurable frequencies, their lifetimes are always less than 15; that is, they survive a maximum of 15 consecutive frames (300 ps). The three main hot spot residues for PAR@ADRP binding are Phe154, Val47, and Ile129, accounting for 20% of the total binding energy ( Figure  6c ). For SIM@ADRP, the electrostatic energy contribution is À23.3 kcal mol À1 owing to the formation of H bonds between the cyclopropylsulfonyl group and Val47 and Ile129 (Figures  3d and 7b) . Although the frequencies of these H bonds are similar to those in the case of PAR@ADRP, their maximum lifetime is longer: 21 (420 ps) for Val47 and 19 (380 ps) for Ile129; this implies a greater electrostatic contribution in the case of SIM@ADRP (Figure 3d ). Ile129 was assigned as the main residue contributing to the SIM@ADRP binding, followed by Leu158 and Pro134 (Figure 6d ).

The phenanthridine group in PAR and the cyclopropylsulfonyl group in SIM are essential for the binding of these drugs to ADRP and Mpro. In the ADRP complexes (Figure 8) , the phenanthridine and thiazole-quinoline groups are oriented in the same direction, and the phenanthridine group is considerably more exposed to the solvent because its large size makes its intercalation into the protein cleft difficult. Figure 8 shows that the methylpyrazine group in PAR occupies the binding site where SIM forms H bonds. The H bonds in SIM@ADRP are not sufficient to overcome the hydrophobic contacts in PAR@ADRP (Figure 3 ), which is more stable according to Table 2 . The MD simulations suggest that although both SIM and PAR interact with ADRP, PAR is a stronger inhibitor of ADRP due to the non-electrostatic contributions.

PAR and SIM weakly interact with Nsp12, with the binding energies of only À16.9 and À18.5 kcal mol À1 , respectively ( Table 2 ). The ligand-receptor contacts are shown in Figure 9 , where the low frequency of contacts between SIM and Nsp12 can be noticed. Only Pro843 and Asn844 have a frequency of contacts >20% in SIM@Nsp12, and Leu789 participates in the formation of H bonds with SIM in 12% of the simulation time. This indicated that, although SIM was docked inside the binding site at the beginning of the 200 ns MD simulation ( Figure  S3a) , it "adsorbed" on the receptor surface via the thiazolequinoline group (Figure 10a and b) .

The PAR@Nsp12 binding mode involves frequent contacts between the cyclopropylsulfonyl group and Tyr516-Arg525; some of these contacts are characterized as H bonds between PAR and Arg523 (5%), Ser519 (8%), and Ala520 (28.5%). In addition to high frequency, these H bonds have relatively long lifetimes: 40 (800 ps) for Ser519 and 33 (660 ps) for Ala520, which reflect their favorable electrostatic contributions. Hydrophobic contacts between the phenanthridine group and Arg525 and Lys591-Asp593 and between the methylpyrazine group and Tyr516, Glu781, and Arg806 strengthen the binding of PAR to Nsp12. Figure S3b shows that, in the case of PAR@Nsp12, the position of PAR does not significantly change along the simulation trajectory, which is contrary to the case of SIM@Nsp12 ( Figure S3a ). For Nsp12 complexes, the effect of solvation energy was considerable. As shown in Table 2 , PAR@Nsp12 is substantially more stable than SIM@Nsp12 in the gas phase (-64.9 versus À39.2 kcal mol À1 , respectively), but in the aqueous phase, SIM@Nsp12 is slightly more stable than PAR@Nsp12. However, the stability order of the four most stable complexes was the same in both phases: SIM@Mpro > PAR@Adrp > PAR@Mpro > SIM@Adrp.

Finally, an additional discussion on the binding energies based on Poisson À Boltzmann/solvent accessibility (PB/SA) and entropic data is provided. Figure 11 depicts the binding energies of the four main complexes evaluated using both the GB/SA and PB/SA approaches.

These approaches differ in the way the solvation energy is obtained (grey bars in Figure 11 ) (Genheden et al., 2015; Wang et al., 2019) . Basically, the solvation energy is divided into two contributions: electrostatic (el) and non-electrostatic (nel); the nel solvation energy is calculated as an empirical linear function of SASA in both models. Figure 10 . Representative snapshots for SIM@Nsp12 (a,b) and PAR@Nsp12 (c,d) complexes. The residues shown were selected according to the relative contact frequency normalized to 100. The arrows indicated the H bonds and the residues involved.

In the GB model, the el solvation energy is analytically calculated using Eq. (2),

where q is the atomic charge, e is the solvent dielectric constant, and f GB is a function of r ij (the distance between the atoms i and j, and the so-called effective Born radii R i and R j Þ: The GB approach is considerably easier to implement and faster than the PB approach, in which the el solvation energy is obtained from the numerical solution of the Poisson equation (3) (Wang et al., 2019) .

where e r ð Þ, U r ð Þ, and q r ð Þ are the dielectric function, electrostatic potential, and solute charge density as a function of the atomic coordinate r, respectively.

As shown in Figure 11 , DG solv is positive and mainly arises from the "desolvation" of the ligand when the solvation energies for the complex and receptor are of the same magnitude and almost cancel out. For the GB/SA approach, DG solv ranges from 32.2 to 39.8 kcal mol À1 , which decreases the stability of the complex by $50% (Figure 11a ). This effect is significantly higher in the case of the PB/SA approach, where DG solv is larger than DG gas (Figure 11b ), making the final binding energy in solution positive. The relative order of complex stability slightly differs between these methods; however, in both cases, the best inhibitor of Mpro is SIM (SIM@Mpro). An attempt to include entropy variation was also made using quasi-harmonic approximation. The decrease in entropy was substantially high for all complexes, with -TDS in the range of 105 À 126 kcal mol À1 (almost twice the value of DG gas ), and the final absolute binding free energy was not considerable. The significant entropy decrease is because of the highly flexible ligands, and, thus, the application of the quasi-harmonic approximation may not be suitable. The applications of the GB(PB)/SA approaches are hindered by many limitations, such as the force-field uncertainties, difficulty in entropy calculation, limited number of MD replicas, and short trajectory length (Wang et al., 2019) . The overall conclusion about all these factors is still unclear, and, therefore, the impacts of these factors should be evaluated and considered during the MD simulations. Herein, the GB/ SA data were chosen, and the relative stability of the complexes was analyzed. When the binding energy was investigated without entropy, the best candidate was SIM@Mpro in the gas and solution phases, regardless of whether the GB or PB approach was used for calculating the solvation energy.

A set of 24 market drugs, including 22 antivirals, were assessed against non-structural proteins (Nsps) of SARS-CoV-2, responsible for the COVID-19 disease. A sequential docking-molecular dynamics (MD) approach was used to screening for the best candidates. The results suggested the paritaprevir (PAR) as a promising antiviral to block most of important functional proteins of SARS-CoV-2. This drug binds very favorable to five of the seven Nsps, including the ADP-ribose-1"-monophosphatase (ADRP), which plays a role in the virus replication process. The complex PAR@ADRP was the second most stable among the six most stable complexes evaluated and the high stability is provided by hydrophobic contacts of the phenanthridine and methylpyrazine rings with the protein binding site. The hot spot residues in ADRP were Phe154, Val46 and Leu124, which are H bonded to the PAR. Besides, Val47 and Ile129 also contribute favorably for PAR@ADRP stability. Overall, the multi-target characteristic of PAR is due to the its structural flexibility, mainly the aromatic side chains, phenanthridine and methylpyrazine, which allow to adjust in the binding site favoring the contacts at the anchor points. Added to the multi-target property of the PAR, the results indicated the simeprevir (SIM) as the best inhibitor of Mpro, a protein active in the viral replication. The high stability of SIM@Mpro comes mainly from non-electrostatic terms, as found for PAR@ADRP, but with important contributions from H bonds with Gly143, Cys145, Asn142 and Ser144 residues, which form a hydrogen bond network around the cyclopropylsulfonyl group of SIM. In summary, the present study allowed to understand the behavior of commercial antivirals in contact with functional proteins of SARS-CoV-2, and suggests that a combination of PAR and SIM could lead to a cooperative biological response against the SARS-CoV-2 virus.

The Supporting Information provides the docking parameters used for virtual screening (Table S1) , docking scores for all drugs evaluated against the seven Nsps (Table S2) , docking scores for the best ten ligands (Table S3) , and the numbering schemes used for all proteins as discussed in the text (Table  S4 ). The structures for all compounds analyzed ( Figure S1 ), comparison of GB/SA binding energies calculated at 60 ns and 200 ns of MD simulation ( Figure S2 ), and the structures of PAR and SIM complexes with Nsp12, showing the large amplitude movement of the SIM ligand along the MD trajectory ( Figure S3 ).

Pharmacoinformatics and molecular dynamics simulation studies reveal potential covalent and FDA-approved inhibitors of SARS-CoV-2 main protease 3CL pro

COVID-2019: The role of the nsp2 and nsp3 in its pathogenesis

Drug interactions with the direct-acting antiviral combination of ombitasvir and paritaprevir-ritonavir

Molecular dynamics with coupling to an external bath

Modeling infectious disease dynamics

Chloroquine and hydroxychloroquine as available weapons to fight COVID-19

Ribavirin, Remdesivir, Sofosbuvir, Galidesivir, and Tenofovir against SARS-CoV-2 RNA dependent RNA polymerase (RdRp): A molecular docking study

Comparative protein structure modeling using MODELLER

Repositioning of 8565 existing drugs for COVID-19

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

A SARS-CoV-2-human protein-protein interaction map reveals drug targets and potential drug-repurposing

Accelerating drug development through repurposed FDA-approved drugs for COVID-19: Speed is important, not haste

Repurposing of known anti-virals as potential inhibitors for SARS-CoV-2 main protease using molecular docking analysis

VMD: Visual molecular dynamics

Triple combination of interferon beta-1b, lopinavir-ritonavir, and ribavirin in the treatment of patients admitted to hospital with COVID-19: An open-label, randomised, phase 2 trial

In silico exploration of the molecular mechanism of clinically oriented drugs for possibly inhibiting SARS-CoV-2's Main Protease

Comparison of simple potential functions for simulating liquid water

Targeting SARS-CoV-2: A systematic drug repurposing approach to identify promising inhibitors against 3C-like proteinase and 2 0 -O-ribose methyltransferase

Crystal structure of Nsp15 endoribonuclease NENDOU from SARS-COV -2

Processing of the SARS-CoV pp1a/ab nsp7-10 region

Ombitasvir/ paritaprevir/ritonavir þ dasabuvir and ribavirin associated druginduced liver injury and syndrome of inappropriate secretion of antidiuretic hormone: A case report

Discovery of new hydroxyethylamine analogs against 3CLpro protein target of SARS-CoV-2: molecular docking, molecular dynamics simulation, and structure-activity relationship studies

Potent antiviral activities of the direct-acting antivirals ABT-493 and ABT-530 with three-day monotherapy for Hepatitis C virus genotype 1 infection

Nsp3 of coronaviruses: Structures and functions of a large multi-domain protein

Respiratory virus shedding in exhaled breath and efficacy of face masks

Crystal Structure of the SARS-CoV-2 Non-structural Protein 9

Research and development on therapeutic agents and vaccines for COVID-19 and related human coronavirus diseases

Genomic characterisation and epidemiology of 2019 novel coronavirus: Implications for virus origins and receptor binding

Drug repurposing using computational methods to identify therapeutic options for COVID-19

ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB

MMPBSA.py : An efficient program for end-state free energy calculations

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

Peptide-like and small-molecule inhibitors against Covid-19

Molecular docking and dynamic simulations for antiviral compounds against SARS-CoV-2: A computational study

Structural and biochemical characterization of the nsp12-nsp7-nsp8 core polymerase complex from SARS-CoV-2

PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data

Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes

silico studies on therapeutic agents for COVID-19: Drug repurposing approach

Glecaprevir and Maraviroc are high-affinity inhibitors of SARS-CoV-2 main protease: Possible implication in COVID-19 therapy

One severe acute respiratory syndrome coronavirus protein complex integrates processive RNA polymerase and exonuclease activities

Targeting SARS-COV-2 non-structural protein 16: A virtual drug repurposing study

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

Synchronization of trajectories in canonical molecular-dynamics simulations: Observation, explanation, and exploitation

Estimates of the severity of coronavirus disease 2019: A model-based analysis. The Lancet Infectious Diseases

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

End-point binding free energy calculation with MM/PBSA and MM/GBSA: Strategies and applications in drug design

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

Analysis of therapeutic targets for SARS-CoV-2 and discovery of potential drugs by computational methods

Crystal Structures of Two Coronavirus ADP-Ribose-1 00 -monophosphatases and their complexes with ADP-Ribose: A systematic structural analysis of the viral ADRP Domain

Real-life results of treatment with ombitasvir, paritaprevir, dasabuvir, and ritonavir combination in patients with chronic renal failure infected with HCV in Turkey

Structural basis for inhibition of the RNA-dependent RNA polymerase from SARS-CoV-2 by remdesivir

The proteins of severe acute respiratory syndrome coronavirus-2 (SARS CoV-2 or n-COV19), the Cause of COVID-19

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

We would like to express our gratitude to the financial support from the Brazilian Agencies CNPq, FAPEMIG and CAPES that provided continuous support to our laboratories. This is a project conducted by members of the Rede Mineira de Qu ımica (RQ-MG).

No potential conflict of interest was reported by the authors.

Fundação de Amparo a Pesquisa do Estado de Minas Gerais.