Remdesivir and Ledipasvir among the FDA-Approved Antiviral Drugs Have Potential to Inhibit SARS-CoV-2 Replication

Abstract

The rapid spread of the virus, the surge in the number of deaths, and the unavailability of specific SARS-CoV-2 drugs thus far necessitate the identification of drugs with anti-COVID-19 activity. SARS-CoV-2 enters the host cell and assembles a multisubunit RNA-dependent RNA polymerase (RdRp) complex of viral nonstructural proteins that plays a substantial role in the transcription and replication of the viral genome. Therefore, RdRp is among the most suitable targets in RNA viruses. Our aim was to investigate the FDA approved antiviral drugs having potential to inhibit the viral replication. The methodology adopted was virtual screening and docking of FDA-approved antiviral drugs into the RdRp protein. Top hits were selected and subjected to molecular dynamics simulations to understand the dynamics of RdRp in complex with these drugs. The antiviral activity of the drugs against SARS-CoV-2 was assessed in Vero E6 cells. Notably, both remdesivir (half-maximal effective concentration (EC50) 6.6 μM, 50% cytotoxicity concentration (CC₅₀) > 100 µM, selectivity index (SI) = 15) and ledipasvir (EC50 34.6 μM, CC₅₀ > 100 µM, SI > 2.9) exerted antiviral action. This study highlights the use of direct-acting antiviral drugs, alone or in combination, for better treatments of COVID-19.

1. Introduction

Coronaviruses are enveloped viruses belonging to the Coronaviridae family and are notorious for their ability to cause respiratory infections in humans [1]. Coronaviruses are categorized into four subgenera (based on differences in protein sequences) named alpha (α), beta (β), gamma (γ), and delta (δ) [2]. To date, β-coronaviruses have been identified as the causative agents of respiratory infections in humans, including the common cold. In 2019, the newly identified coronavirus SARS-CoV-2 causing COVID-19 caused a global pandemic. The genome of coronaviruses has a size ranging from 27 to 32 kilobase pairs (kbp) and codes for open reading frames 1a and 1b (ORF1a and ORF1b), i.e., polyproteins that regulate viral replication [3,4]. Four other coronavirus strains have been known previously to cause disease in humans and are named 229E, HKU1, NL63, and OC43 [5].

There are different nonstructural proteins (nsps) named nsp1 through nsp16 in coronaviruses that play a major role in genome transcription and replication; however, the exact function of a few nsps remains unidentified. Structural proteins play an important role in viral infection and virion assembly, whereas the spike-like S protein facilitates surface attachment to host cells [5,6]. The M protein with transmembrane domains binds to the nucleocapsid and shapes the virion [7,8]. The envelope protein (E protein) is important for viral infection pathogenesis because it performs a crucial function in virion budding and assembly. The N protein is composed of two domains: one with the ability to bind to the viral genome and a nsp3 triggering the replicase–transcriptase complex and the encapsulation of the viral genome. Most RNA viruses, with the exception of retroviruses, require an RdRp for transcription and replication of their genome; therefore, this protein is crucial for their survival [9]. RdRp is highly conserved among RNA viruses and is composed of a 900-amino-acid catalytic core that consists of the palm, fingers, and thumb domains [10]. RdRp has been identified as a conserved protein within a group of RNA viruses and therefore could be an attractive research object for analyzing the process of nucleic acid synthesis and for the development of antiviral drugs [10,11].

Soon after the SARS-CoV-2 pandemic, many in vitro studies and clinical trials were conducted to identify potent drugs with high antiviral activity. In this regard, favipiravir is one of the first antiviral drugs authorized by the National Medical Products Administration of China for use against SARS-CoV-2 [12]. In the current scenario, the efficacy of the antiviral drugs may not be as effective against virus due to single nucleotide polymorphisms which changes the viral proteins [3]. A similar phenomenon can also be observed in SARS-CoV2, whereby the viral proteins (i.e., spike protein) acquire mutations due to SNPs and consequently escape being targeted by antiviral drugs [13,14,15,16]. However, the conserved protein such as RdRp could be an appropriate target to understand the mechanism of nucleic acid synthesis and the subsequent identification of already available antiviral drugs [17,18,19,20,21]. Several other drugs, such as remdesivir, ribavirin, and sofosbuvir, are currently in different phases of clinical trials for evaluation of their anti-RdRp activity. Nonetheless, a multistep approach is needed to speed up the drug discovery process [12,22,23,24].

In this context, to overcome the problems associated with the discovery of effective antiviral drugs, the current study was conducted to determine potential effectiveness of FDA-approved antiviral drugs against SARS-CoV-2 replication. Here, we employed a combined multistep in vitro and in silico approach, involving virtual screening, docking, and molecular dynamics (MD) simulations to identify a lead compound effective against the virus. The best drug was then selected based on a binding score, stability, and its interaction with key amino acid residues to be confirmed by in vitro tests. The antiviral activity of the top hits (remdesivir, ledipasvir, and paritaprevir) was then assessed in Vero E6 cells. The high antiviral activity of remdesivir (half-maximal effective concentration (EC₅₀) = 6.6 μM) and ledipasvir (EC₅₀ = 34.6 μM) observed in Vero E6 cells against SARS-CoV-2 highlights their potential utility as an effective therapeutic intervention to be tested in combination to achieve better efficacy for further clinical evaluation. The findings of this study will also help physician scientists to treat COVID-19 patients with an appropriate drug.

2. Materials and Methods

2.1. Protein Model Preparation and Active-Site Identification

The three-dimensional (3D) structure of the SARS-CoV-2 RdRp protein was retrieved from the Protein Data Bank (PDB; ID: 6M71) [25]. The retrieved cryo-electron microscopy structure was prepared using the Molecular Operating Environment (MOE) software [26]. All water molecules and associated cofactors (nsp7–8) were removed, and hydrogen atoms were added to the structure. All missing atoms were modeled, and energy was minimized to remove any steric clashes and to optimize bond lengths and angles. For energy minimization, the AMBER10: EHT force field was utilized with a 0.01 rms kcal/mol gradient. The structural topology was reviewed for any abnormalities. The active site of RdRp is located in the seven conserved motifs (A to G). The sterol-sensing domain sequence (amino acid residues 759–761, particularly residues K545 and R555) was found to be a hotspot.

2.2. Library Preparation

All the available antiviral drugs (n = 63) approved by the FDA in the last 50 years were selected for virtual screening, and their 3D structures were downloaded from the PubChem database [27,28]. The antiviral drug library was prepared in MOE, explicit hydrogen atoms were added, and all molecules were subjected to energy minimization by means of the Merck molecular force field 94× (MMFF94×) with a root mean square gradient of 0.1 [29,30].

2.3. Virtual Screening and Molecular Docking

Docking was performed using Moe-Dock (Chemical Computing Group Inc., Montreal, QC, Canada). Each FDA-approved antiviral drug was docked separately to the active site of the RdRp protein. At the first step, water molecules were removed, and 3D protonation and energy minimization were performed using the MOE software with the following parameters: force field MMFF94X, gradient 0.05, and current geometry. The FDA-approved antiviral drug library was docked into the catalytic binding site of RdRp (Arg555, Val557, Asp618, Asp623, Thr680, Asn691, Ser759, Asp760, and Asp761). Next, active-site residues were selected, with 10 conformations, and the prepared library was screened against the selected residues. The protein and ligand atoms were kept flexible, and the triangle-matcher placement method was used along with the London dG scoring function. Among the top 10 hits, three drugs were selected (based on a docking score, binding mode, hydrogen bonding, and interactions) for MD simulations; namely, remdesivir, ledipasvir, and paritaprevir were chosen, and their intermolecular interactions were studied by means of PyMOL and MOE [31]. After that, the docked complexes were subjected to MD simulations.

2.4. The MD Simulation Protocol

2.4.1. Ligand Topology Generation

The official CHARMM general force field (CGenFF) server was used to build the topology of the selected ligands for the MD simulations [32]. Initially, the ligand molecules were converted into mol2 format, hydrogen atoms were added into the correct tautomeric and protonation state, and the bond order was corrected accordingly. Next, the ligands in the mol2 format were uploaded to the CGenFF server to generate a CHARMM-compatible stream file. The stream file contains ligand topological information, such as atom types, charges, and bond parameters. The stream file was then converted into a Gromacs-compatible file using a python script that generates a PDB file, parameters, and topology files for the ligand. The generated ligand molecular structure was then saved in the Gromacs-compatible format.

2.4.2. MD Simulations

The RdRp complexes with the chosen FDA-approved drugs were subjected to MD simulations to discern the dynamics and interaction behavior of these compounds. The all-atom simulation method was used to gain insights by solving Newton’s equation of motion. Initially, SARS-CoV-2 RdRp in complex with a drug (remdesivir, ledipasvir, and paritaprevir) and the apo-RdRp protein (PDB ID: 7BV1) were subjected to energy minimization to remove steric effects and to optimize the structure. The CHARMM36 all-atom force field was utilized [33]. Periodic boundary conditions were applied, and the protein was placed inside a cubic box with 10 Å distance between the boundary and the protein surface. The TIP3P water model was employed to solvate the systems, and the system was neutralized by the addition of Na⁺ ions. The preprocessing steps were then executed, which included energy minimization and temperature and pressure equilibration. Energy minimization was carried out by the steepest-descent algorithm for the subsequent 50,000 steps, followed by NVT and NPT ensemble for 100 ps in order to attain a point of maximum force of 1000 kJ mol⁻¹ nm⁻¹. Then, the temperature equilibration was performed via the V-rescale algorithm at 300 k with a time constant of 0.1 [34] The particle mesh Ewald algorithm [35] was used to calculate long-range electrostatic interactions, while short-range electrostatic and van der Waals were calculated by specifying 1.2 nm cutoff distance. Additionally, the Parrinello–Rahman algorithm was applied to achieve pressure equilibration at 1 bar. After that, 100 ns (nsteps 50,000,000) MD simulations were carried out for protein–ligand complexes and the apo-RdRp protein using the Gromacs software [36]. The data analyses were performed by means of Gromacs built-in tools, MOE, and PyMOL.

2.5. Post-MD Simulation Data Analysis and Visualization

The trajectories obtained throughout the MD simulations were researched by further analyses, such as root mean square deviation (RMSD) computation to measure the system stability, whereas root mean square fluctuation (RMSF) was employed to evaluate flexibility at the amino acid residue level. The radius of gyration (R_g) was calculated to measure structure compactness. Principal component analysis (PCA) was performed to evaluate functional dynamics of the protein. Additionally, to characterize stable and variant states of the protein, free energy landscape (FEL) analysis was performed. Lastly, to identify correlated motions of residues, a dynamic cross-correlation matrix (DCCM) analysis was performed. For all these analyses, Grace (plasma-gate.weizmann.ac.il/Grace), R packages, and Gromacs built-in and standalone tools were used.

2.6. Principal Component Analysis (PCA)

This is a dimensionality reduction method predominantly serving to demonstrate slow and functional motions of biological molecules [37]. To perform PCA, a covariance matrix was calculated by diagonalizing and solving eigenvalues and eigenvectors. The direction of the motion is represented by eigenvectors, whereas the magnitude of the motion and direction is denoted by eigenvalues. Furthermore, the covariance matrix for the illustration of PCA was computed with Gromacs analysis tools g_covar and g_anaeig [38]. Consequently, the trajectories were changed to DCD format by means of stand-alone software Wordom [39]. The PCA was further executed via the Bio3D analysis tool by a method described previously [40].

2.7. The Free Energy Landscape (FEL)

An FEL is determined to characterize all possible conformational changes of a protein in MD simulations [41,42]. The FEL represents two variables that reflect a stable and transient state of the protein and was computed here from a probability distribution composed of the first two eigenvectors from the essential plane. Positions of the interacting molecules in the system were characterized by focusing on their respective energy levels [43]. In the calculation of the FEL, protein stability was identified by Gibb’s free energy calculation. Here, we utilized gmx sham, a Gromacs tool, to compute the FEL. In this study, the apo form and protein–drug complexes were evaluated via the following equation:

G_i = −K_BTln(N_i/N_max)

where K_B is the Boltzmann constant, T denotes temperature (300 K), N_i is the population of bin i, and N_max is the population of the most populated bin. A color-coded model depicts various energy levels.

2.8. The Dynamic Cross-Correlation Matrix (DCCM)

The DCCM analysis was carried out to identify correlated motions of residues upon drug binding, and the Bio3D package available in the R software was employed to calculate residue–residue dynamic cross-correlations [44]. To build the matrix, only Cα atoms were selected. Covariations of the matrices were calculated on calling “cov2dccm” upon calculating Pearson’s covariance matrix correlation coefficients from the coordinates. Based on the following equation, the cross-correlation ratio and matrix (C_ij) represent time-correlated data between atoms i and j of a protein [45,46]:

$${\mathrm{C}}_{\mathrm{ij}}={⟨\mathsf{\Delta }{\mathrm{r}}_{\mathrm{i}}\cdot \mathsf{\Delta }{\mathrm{r}}_{\mathrm{j}}⟩∕\{⟨\mathsf{\Delta }{\mathrm{r}}_{\mathrm{i}}⟩2⟨\mathsf{\Delta }{\mathrm{r}}_{\mathrm{i}}⟩2\}}^{1/2}$$

In the above equation, both Δr_i and Δr_j show average locations of the ith and jth residues, respectively. The angular brackets represent the average time. The positive values indicate correlated motions, and negative values indicate atomic displacements in opposite directions.

2.9. MM-PBSA Calculation

The MM-PBSA method was used to calculate the binding free energy between protein and ligand complexes. In total, 100 frames were extracted from the trajectories and the total energy of the system was calculated through the following molecular Mechanic/Poisson–Boltzmann Surface Area (MM-PBSA) equation

$$\mathsf{\Delta }\mathrm{G}\left(\mathrm{Binding}\right)= \mathrm{G}\left(\mathrm{Complex}\right)- \mathrm{G}\left(\mathrm{Receptor}\right)- \mathrm{G}\left(\mathrm{Ligand}\right)$$

In the above equation, G (complex) represents the total free energy of the protein-ligand complex, whereas G (receptor) and G (ligand) are total free energies of the isolated ligand and protein in solvent. Thus, the ‘g_mmpbsa’ tool from the Gromacs package was used to perform MM-PBSA calculations [47]. A detailed explanation of the method used in the calculations can be found in previous studies [48,49,50]

2.10. In Vitro Activity of the Antiviral Drugs

The compounds were purchased from MolPort (https://www.molport.com), and 20 µM stock solutions were prepared in dimethyl sulfoxide (50–100 µL). The sequence of a SARS-CoV-2 strain isolated from COVID-19 patients by Korea Disease Control and Prevention Agency was used (BetaCoV/Korea/KCDC03/2020). One-hour post infection, the viral inoculum was removed, and Vero E6 cells were treated with serial dilutions of the candidate drugs in the infection medium. After 48 h, the cells were stained with an antibody against the constituent protein of the virus to evaluate the degree of inhibition of viral replication in the sample, and cell viability (drug toxicity) was assessed by examining the presence of the nucleus in the cells. If a drug was effective in a sample, all the cell nuclei were stained (drug toxicity 0%), and the virus component protein was not stained by the antibody (inhibition 100%). Dilutions of the solvent in the infection medium were utilized to set up mock-treated controls. The virus-inhibitory effect was checked using the antibody against the viral protein as a vehicle, and cell viability (drug toxicity) was evaluated via the nuclear staining of cells.

3. Results and Discussion

3.1. Virtual Screening and Molecular Docking

In this study, we chose the virtual-screening approach to identify potential antiviral drug candidates against SARS-CoV-2 RdRp that target active-site residues (Figure 1). The active site of RdRp is highly conserved among various microbes [51]. The criteria for selecting the best compounds are the docking score (S) and multiple interactions with the active-site residues. The top 10 compounds that satisfied the filtering criteria were visualized manually. Finally, out of these top hits, three antiviral drugs were selected for comparative analysis.

3.2. The Analysis of the Interaction of Top Hits with RdRp

The molecular docking of FDA-approved antiviral drugs with the SARS-CoV-2 RdRp crystal structure was performed next. The docking scores of the top hits (remdesivir, ledipasvir, and paritaprevir) from the drug library implied a strong interaction with the key residues. For ledipasvir, the docking score was found to be −8.9 kcal/mol. It forms three hydrogen bonds with active-site residues. The amidine group of R555 engages in two hydrogen-bonding interactions with the imidazole group of ledipasvir, with bond lengths of 2.1 and 2.6 Å. Similarly, the carboxyl group of another amino acid, V166, forms a hydrogen bond with the methoxycarbonyl amino group (bond length of 2.2 Å;

Table 1. Docking scores and hydrogen-bonding interactions with important binding-site residues.

Drug Name	Hydrogen Bonds	Bond Length (Å)	Force Field	MOE Docking Score (kcal/mol)
Ledipasvir	R555	2.1 and 2.6	CHARMM36	−8.9
	V166	2.2
Remdesivir	T556	2.1 and 2.4	CHARMM36	−8.3
	D760	2.3 and 2.5
	N691	2.5
Paritaprevir	D761	2.1	CHARMM36	−9.4
	S814	2.0

, Figure 2a).

Remdesivir, with a docking score of −8.3 kcal/mol, yielded five hydrogen bonds. The carboxyl group of T556 forms two hydrogen bonds with the oxygen atom of the 5-cyano-3,4-dihydroxyoxolan-2-yl moiety, with bond lengths of 2.1 and 2.4 Å. Similarly, the deprotonated carbocyclic group of D760 and amino group of N691 form two and one hydrogen bonds, respectively, with the hydrogen and nitrogen atoms of the nucleobase substituent 4-aminopyrrolo[2,1-f][1,2,4]triazin-7-yl group, with bond lengths of 2.3 and 2.5 Å (Table 1, Figure 2b).

The docking score of paritaprevir was −9.4 kcal/mol, and this drug participates in hydrogen bonding with the key binding-site residues (D761 and S814). The carboxyl group of D761 forms a hydrogen bond with the amide group of an aza-macrocyclic complex with a bond length of 2.1 Å. Additionally, the amino group of S814 engages in a hydrogen bonding with the carbonyl group attached to the acyl-sulfonamide moiety with a bond length of 2.0 Å (Table 1, Figure 2c). This study suggests that remdesivir engages in stronger interactions with RdRp than ledipasvir and paritaprevir do, thereby inhibiting the template entry site more effectively. It is known that aspartate, serine, and arginine play critical roles in the formation of hydrogen bonds [52]. Furthermore, the number of hydrogen bonds between the protein and ligands is shown in (Figure 2a–c). During MD simulations, the number of formed hydrogen bonds in both RdRp–ledipasvir and RdRp–remdesivir complexes varied between 0 and 6, while for RdRp–paritaprevir, the hydrogen bond number varied between 0 and 7.

3.3. Measurement of Variations in the Cα Atoms of RdRp in the Presence and Absence of Drugs

All the RdRp–ligand complexes were simulated in an aqueous environment for 100 ns each to measure the conformational changes and discern the dynamics of stability. The deviation of backbone atoms was measured by means of RMSD. RMSDs of all the systems are given in Figure 3. To characterize the dynamic behavior of RdRp during the interaction with the selected FDA-approved antiviral drugs (remdesivir, ledipasvir, and paritaprevir), MD simulations of the RdRp-drug complexes, with each selected drug, were performed in an explicit water environment. Dynamic behaviors of the protein and ligands were analyzed individually using MD simulation trajectories. The RMSDs of the complexes were assessed and compared with the RMSD of apo-RdRp (Figure 3a). In this analysis, the RdRp–ledipasvir complex displayed steady incremental deviation during the simulation with RMSD starting from 2.1 Å and reaching 3.3 Å. Similarly, RdRp–remdesivir showed incremental deviation (2 to 3.6 Å) throughout the simulation, along with some acceptable convergence during the time intervals. By contrast, the RdRp–paritaprevir complex manifested a dramatic increase in the deviation from 2 to 4 Å between nanosecond 0 and nanosecond 6, but soon the complex stabilized and stayed stable throughout the rest of the simulation, with the RMSD reaching 4.3 Å. In the comparison, the RMSD of RdRp in the complexes with remdesivir and ledipasvir yielded a trajectory similar to that of the apo-form of RdRp, with slight variation during the time interval. Overall, the RdRp–ligand complexes attained an average RMSD of approximately 3.7 Å by the end of the simulations.

Moreover, we measured residue flexibility by means of RMSF. In the case of apo-RdRp, greater fluctuation primarily occurs in residues K50 (3.099 Å), R116 (7.783 Å), N403 (4.14 Å), and D910 (6.617 Å), which are present in loop regions, while no significant fluctuations were observed in other regions (Figure 3b). Additionally, RMSFs of the active-site residues, R555 (0.876 Å), V166 (1.301 Å), T556 (0.854 Å), D760 (1.157 Å), N691 (0.727 Å), D761 (1.038 Å), and S814 (1.01 Å), indicated no significant fluctuations (Figure 3b).

Nevertheless, for the RdRp–remdesivir complex, significant fluctuations were noted in the loop regions, especially surrounding residues Y69 (6.38 Å), F103 (5.317 Å), and T262 (9.17 Å). On the other hand, the scores (Cα RMSFs) for the catalytic binding-site residues T556, D760, and N691 were 1.512, 1.536, and 0.94 Å, respectively. Additionally, a similar fluctuation pattern was identified in the loop region of the RdRp–ledipasvir complex; in particular, residues Y69, F102, H362, M906, and E919 featured greater fluctuations. The scores (Cα RMSFs) for the active-site residues R555 and V166 were 1.741 and 1.34 Å, respectively.

In the RdRp–paritaprevir complex, loop region residues Y69 (5.307 Å), F102 (4.221 Å), T262 (3.157 Å), R365 (3.156 Å), E431 (3.806 Å), K508 (3.67 Å), D824 (3.992 Å), G852 (4.02 Å), and L907 (5.801 Å) showed greater fluctuation. The scores (Cα RMSFs) for the catalytic binding-site residues D761 and S814 were 1.489 and 1.791 Å, respectively.

Overall, the protein–ligand interactions were found to differently affect residue flexibility and internal dynamics. The average RMSF for the RdRp–drug complexes are 1.49 Å, which is slightly higher than that of unbound RdRp (1.00 Å).

3.4. Calculation of Equilibrium Conformation of Systems

Structural compactness of all the systems was measured through the calculation of R_g of their MD trajectories, and average values were obtained. The effects of drug binding on R_g of the protein were measured and compared with the apo-RdRp protein structure. The average R_g of apo-RdRp turned out to be 29.9 Å, and the system remained stable and compact. For the RdRp–ledipasvir complex, the average R_g was 30.2 Å and showed negligible deviation in terms of compactness and stability as compared to the apo-structure. Similarly, an average R_g of 30.0 Å was observed in the RdRp–paritaprevir complex. On the contrary, for the RdRp–remdesivir complex, R_g increased to 30.5 Å (Figure 3c). Thus, RdRp in all the complexes showed sustained stability and compactness corresponding to native RdRp.

3.5. Analysis of Essential Dynamics of Protein

Essential dynamics in a protein are controlled by switching between different conformations, and the phenomena governing this modular nature of the protein are controlled by overall collective motions [53]. This concept is important in many biological processes and plays an important role in biological signaling pathways. In this context, considerable flexibility and rigidity are required for a protein to be functional, especially in a binding site [54,55]. Moreover, a strong interaction may restrict protein movement, thereby affecting its biological activity [56]. Therefore, PCA was carried out to investigate the collective motion of unbound and bound drugs (ledipasvir, remdesivir, and paritaprevir) in the MD trajectories. In this dimensionality reduction method, a projection of two principal components, PC1 and PC2, is calculated by diagonalizing the covariance matrix of eigenvectors to describe the subspace in which maximum protein dynamics occur. Furthermore, by PCA, the MD trajectories of all the protein–ligand complexes and of the apo-protein were examined to evaluate the conformational and structural changes upon ligand binding. The dynamic motions of the protein determined by PCA are depicted in Figure 4. The unbound protein (apo-RdRp) clearly possesses fewer stable clusters compared to the RdRp–drug complexes, particularly in PC2. In the RdRp–drug complexes, a greater number of stable clusters in some protein regions were identified, indicating a reduction in the overall collective motion because the drug binding effectively limits the protein backbone movements in certain regions. The conformational space covered in the case of the ledipasvir and remdesivir complexes proved to be broader than that of the complex with paritaprevir. These results mean that upon ligand binding, the overall conformation of the RdRp protein changes, which may alter the desired protein functionality.

3.6. Protein Folding Dynamics Exploration

The landscapes of energy minima of both apo-RdRp and drug-bound RdRp were visualized using the FEL against two principal components, PC1 (RMSD) and PC2 (R_g), which showed ΔG values between 0 and 10 kJ/mol (Figure 5). The stability of the protein is represented by the darkest and centralized blue regions of the FEL. This region in the plot denotes the energy minima of different conformations and represents the stability of the protein complex. As illustrated in Figure 5a–d, the lowest-energy state of the FEL of apo-RdRp is achieved after 13 ns, whereas the RdRp in complex with remdesivir, ledipasvir, or paritaprevir achieves stability at 15, 77, and 35 ns, respectively. Moreover, the low-energy zones of the apo form and of RdRp in complex with remdesivir or ledipasvir are larger than those of the paritaprevir complex. This finding suggests that the RdRp–paritaprevir complex goes through a comparatively lengthy transition state to achieve an equilibrium. Nonetheless, overall, the selected drugs (remdesivir, ledipasvir, and paritaprevir) have a potential to cause the RdRp protein to enter a local energy minimum state.

3.7. Time-Correlated Protein Domain Motions and Molecular Flexibility

The DCCM is a 3D matrix that graphically depicts time-correlated information on correlations among a protein’s amino acid residues. The most common method used to analyze residues based on time-correlated data is visual pattern recognition. In this regard, DCCM analysis of Cα atoms was performed throughout the simulations to probe the dynamics of apo-RdRp and of RdRp in complex with a drug (remdesivir, ledipasvir, or paritaprevir), as illustrated in Figure 6. The residues with highly positively correlated motions are shown as red regions, while highly inversely correlated motions are blue regions. According to Figure 6a, the active-site residues in the apo form undergo stronger correlated motions. In the RdRp–remdesivir complex, the active-site interacting residues (T556, N691, and D760) showed a positive correlation, just as in the apo-form. Similarly, interacting residues V166 and R555 in the ledipasvir complex and D761 and S814 in the paritaprevir complex also manifested strong correlations. These results uncovered the presence of stronger cross-correlation dynamics between residues in the RdRp–drug complexes (Figure 6b–d), suggesting strong interactions and better stability. As documented in this study, residues in these regions also have similar RMSF values in both the apo form and RdRp bound with the selected drugs. Thus, it can be concluded that the pattern of correlation observed here improves overall system stability.

3.8. Binding Free Energy Calculations

The MM-PBSA method was employed to calculate the binding free energy from the trajectory obtained during MD simulation. For each protein–ligand complex, binding free energy (∆G_bind), van der Waals energy, electrostatic energy, and polar solvation energy were calculated as shown in Figure 7a. It has been observed in previous studies that binding free energy values lower than −30 kJ/mol can be taken for binding, but lower binding free energy values are considered to be more favorable for the interaction [57,58]. In the current study, RdRp-drugs (remdesivir, ledipasvir, and paritaprevir) showed acceptable binding free energy values as shown in Figure 7b. The cumulative binding energy contributed for selected drugs remdesivir, ledipasvir, and paritaprevir were −94.877, −81.18, and −132.108 kJ/mol respectively. In all three complexes, the van der Waals interactions significantly contributes in the binding energy. Overall, our analysis established that the selected drugs have the potential to bind tightly to the SARS-CoV2 RdRp protein.

3.9. In Vitro Antiviral Activity

To determine the effectiveness of the selected FDA-approved drugs against the replication of highly pathogenic SARS-CoV-2, the antiviral activities of the three drugs were evaluated next. In this assay, (a) remdesivir, an adenosine analogue that has been evaluated before the 2018 Kivu Ebola epidemic [59]; (b) ledipasvir, an orally administered NS5A inhibitor used against hepatitis C virus (HCV) [60]; and (c) paritaprevir, a 3/4A protease inhibitor prescribed against HCV infection [61] were evaluated against SARS-CoV-2. Here, we performed antiviral assays on the Vero E6 cell line, and remdesivir served as a positive control according to the observed antiviral activity. We assessed virus replication in the cultured cells exposed to various drug doses (100, 33, 11, 3.7, and 1.2 μM). Judging by the results, 50% cytotoxicity concentration (CC₅₀) of remdesivir was >100 μM, EC₅₀ was 6.6 μM, and the selectivity index (SI) was >15 (Figure 8a). By contrast, ledipasvir manifested antiviral activity with an EC₅₀ of 34.6 μM, and CC₅₀ was >100 μM with SI > 2.9 (Figure 8b). Paritaprevir was found to have an antiviral activity with EC₅₀ of 33.9 μM, and CC₅₀ was 28.5 μM with SI = 0.84 (Figure 8c). Although remdesivir showed a higher effective drug concentration than ledipasvir did, its SI remains of interest.

Taken together, these results mean that remdesivir has a more potent in vitro antiviral activity than ledipasvir does at various doses for the control of viral replication within 24–48 h in our assay system. We can hypothesize that the mechanism is the inhibition of the replicase complex of the virus. Unfortunately, paritaprevir does not have any useful antiviral activity against SARS-CoV-2. Remdesivir is effective against SARS-CoV-2 but has adverse effects [62]. Furthermore, it has a half-life of 0.4 h in nonhuman primates, and for this reason, it offers human angiotensin-converting enzyme 2 (ACE2) inhibition of shorter duration. By contrast, the nucleoside triphosphate metabolite of remdesivir has a half-life of 14 h in nonhuman primates and approximately 20 h in humans [63]. On the other hand, ledipasvir is not extensively metabolized, and its median terminal half-life is 47 h [64]. Due to its longer half-life, ledipasvir shows promise for prolonged action on extracellular human ACE2 as a target. Similarly, the intracellular concentration effective against RdRp in host cells is apparent from a mechanistic insight into known antiviral activity against HCV. Additionally, ledipasvir is a highly protein-bound drug (>99.8%), implying better free-drug concentration at the target site, and therefore should be suitable for an effective therapeutic intervention. Consequently, a multitarget treatment strategy involving a synergistic combination of drugs with different mechanisms of action may offer more robust practical treatment of COVID-19.

4. Conclusions

The present study identified suitable FDA-approved drugs for the inhibition of the SARS-CoV-2 RdRp enzyme. In this regard, a structure-based virtual screening of FDA-approved antiviral drugs was implemented to find lead molecules on the basis of the interaction with the RdRp active site. MD simulations were performed for 100 ns to investigate the dynamic behavior of protein–ligand interactions, revealing that the protein–ligand complex maintains a stable conformation with lower protein–ligand interaction energy. The top hits were then confirmed by an in vitro assay to validate the efficacy of the selected drugs experimentally against SARS-CoV-2 replication. In this context, remdesivir and ledipasvir exerted antiviral action against the virus, and unexpectedly, the antiviral activity of our selected drugs was recently identified in another cell-based screening assay against SARS-CoV-2 [65,66]. This evidence may guide the selection of downstream experiments in further studies and may be a starting point for further confirmation of the selected drugs in a synergistic combination for use in a biologically relevant and more complex preclinical model of COVID-19.