Combined use of the hepatitis C drugs and amentoflavone could interfere with binding of the spike glycoprotein of SARS-CoV-2 to ACE2: the results of a molecular simulation study.

The worldwide rapid spread of the COVID-19 disease necessitates the search for fast and effective treatments. The repurposing of existing drugs seems to be the best solution in this situation. In this study, the molecular docking method was used to test 248 drugs against the receptor-binding domain (RBD) of spike glycoprotein of SARS-CoV-2, which is responsible for viral entry into the host cell. Among the top-ranked ligands are drugs that are used for hepatitis C virus (HCV) treatments (paritaprevir, ledipasvir, simeprevir) and a natural biflavonoid amentoflavone. The binding sites of the HCV drugs and amentoflavone are different. Therefore, the ternary complexes of the HCV drug, amentoflavone, and RBD can be created. For the 5 top-ranked ligands, the validating molecular dynamics simulations of binary and ternary complexes with RBD were performed. According to the MMPBSA-binding free energies, the HCV drugs ledipasvir and paritaprevir (in a neutral form) are the most efficient binders of the RBD when used in combination with amentoflavone.Communicated by Ramaswamy H. Sarma.

Introduction

Beginning in December 2019, the COVID-19 disease, which is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has been rapidly spreading around the world. As of December 26th, there have been more than 79,860,000 cases and 1,750,900 deaths due to COVID-19 (https://coronavirus.jhu.edu/map.html). The phylogenetic analysis of SARS-CoV-2 revealed that it belongs to the Betacoronavirus genus and is closely related to the bat coronaviruses and human SARS-CoV (Lu et al., 2020; Xu et al., 2020). The genome of SARS-CoV-2 is a positive-sense single-stranded RNA. It encodes four structural proteins (spike, envelope, membrane, and nucleocapsid), 16 nonstructural proteins (NSPs), and 9 accessory proteins (Gordon et al., 2020). The NSPs are expressed as two large polyproteins that are proteolytically cleaved into 16 smaller proteins (Arya et al., 2021). To enter the host cell, SARS-CoV-2 uses its spike glycoprotein to bind to the angiotensin-converting enzyme 2 (ACE2) receptor (Zhang et al., 2020). This makes the SARS-CoV-2 spike glycoprotein an attractive target for drug development; specifically, ligands that bind to the spike glycoprotein–ACE2 contact surface could interfere with the ACE2-spike interaction and prevent the penetration of the virus into the cell. When considering that the development of new drugs for COVID-19 could take years, the repurposing of already existing drugs appears to be a good alternative.

Another possible drug targets of SARS-CoV-2 are the NSPs that play an essential role in the virus replication cycle (Rohaim et al., 2021). Among them are the main protease (NSP5) and the papain-like protease (NSP3), which are responsible for the cleavage of the polyproteins into the mature NSP, the RNA-dependent RNA polymerase (NSP12), the helicase-triphosphatase (NSP13), the 3′-5′-exoribonuclease (NSP14), the RNA uridylate-specific endoribonuclease (NSP15), and N7- and 2′-O-methyltransferases (NSP10 and NSP16) (Arya et al., 2021). Due to the efforts of multiple research groups, the structural information on the SARS-CoV-2 proteins is rapidly increasing (Mariano et al., 2020). Currently, there are more than 1000 entries of SARS-CoV-2 proteins deposited in the protein data bank (PDB) (Berman et al., 2000, http://www.rcsb.org/). Recently, the CoV3D database of high-resolution coronavirus protein structures was created (Gowthaman et al., 2021). This provides opportunities for structure-based drug design for COVID-19.

The modern drug discovery process relies on computer simulation methods (Sliwoski et al., 2014). First, the 3D structure of the biological target is taken from the database (PDB). Then the appropriate library of ligands is chosen, and the virtual screening is performed by using the molecular docking method (Pinzi & Rastelli, 2019). The structure of the target is usually kept rigid during the docking simulation. Sometimes the ensemble of target structures is used for docking. The hit ligands are filtered according to the specified criteria (Lipinski’s rule of five, ADMET properties) (Lipinski et al., 2001; Jia et al., 2020). The docking complexes of top ligands are further refined by molecular dynamics (MD) simulation with explicit solvent (Ganesan et al., 2017). During the MD simulation, both the target and ligand could undergo conformational changes to better fit each other. The intermolecular target–ligand interactions could be tracked in real-time (Naqvi et al., 2018). The MD trajectory of the complex can be also used to estimate the ligand-binding free energy using the end-point methods MMPBSA or MMGBSA (Genheden & Ryde, 2015). The results are known to be strongly dependent on the simulation details (force field, charge models, continuum solvation method, solute dielectric constant, atomic radii set, method of entropy calculation, sampling) and, generally, could not be directly compared to the experimental values (Wang et al., 2019). At the same time, the relative affinities of ligands to the same target are usually well reproduced (Wang et al., 2016). Thus, the docking scoring can be improved using MD post-processing (Liu et al., 2018).

Most of the molecular modeling studies against SARS-CoV-2 use the main protease as a drug target (Aanouz et al., 2020; Alamri et al., 2020; Badavath et al., 2020; Bhardwaj et al., 2020, 2021; Ghosh et al., 2020; Jiménez-Alberto et al., 2020; Khan, Ali, et al., 2020; Kumar et al., 2020; Lokhande et al., 2020; Mishra et al., 2020; Saravanan et al., 2020). Perhaps this is because of the pivotal role of this enzyme in the virus life cycle and the high conservation of its active site (Ullrich & Nitsche, 2020). At the same time, MD simulation studies reveal the high conformational flexibility of the binding pocket of the main protease of SARS-CoV-2 (Bzówka et al., 2020; Grottesi et al., 2020). Thus, the use of a single rigid target for docking could be not appropriate in the case of the main protease.

A large number of drug-repurposing studies are focused on the spike–ACE2 interaction surface (Alazmi & Motwalli, 2020; Awad et al., 2020; de Oliveira et al., 2020; Deganutti et al., 2021; Hakmi et al., 2020; Tao et al., 2021; Trezza et al., 2020; Unni et al., 2020; Wei et al., 2020). There are also studies searching for inhibitors for RNA-dependent RNA polymerase (Elfiky, 2020), papain-like protease (Ibrahim et al., 2020), and RNA endoribonuclease (Khan, Jha, Singh, et al., 2020; Sharma et al., 2021). In many works, a multitarget approach was used (Albohy et al., 2020; Chikhale et al., 2020; Khan, Jha, Amera, et al., 2020; Maffucci & Contini, 2020; Mishra et al., 2021; Patil et al., 2021; Puttaswamy et al., 2020; Rameshkumar et al., 2021; Swargiary et al., 2020). Despite the large number of drug-repurposing studies against SARS-CoV-2, there is little agreement about the identified hit candidates and almost no experimental validation of the results (Dotolo et al., 2021). Therefore, additional molecular modeling studies are needed to obtain the consensus results which can be further verified in clinical trials.

The aim of this study was to find the potential inhibitor for the spike glycoprotein that binds in the ACE2-binding region by computer simulation methods. We used the molecular docking method to test the list of 248 BLDpharm drugs against the spike glycoprotein of SARS-CoV-2. This list includes compounds with antiviral activity. Among them are FDA-approved drugs and natural substances used in traditional medicine. All these compounds can be freely purchasable. For the 5 top-ranked drugs, the 120 ns MD simulations of complexes with explicit solvent were performed and MMPBSA-binding free energies were calculated. The binding sites of ligands and the interactions with the spike glycoprotein were analyzed in detail.

Computational methods

Preparation of the receptor and ligand structures

The receptor-binding domain (RBD) of the spike glycoprotein of SARS-CoV-2 was chosen as a docking target. Its coordinates were taken from the cryo-EM structure of the RBD/ACE2-B0AT1 complex deposited in the PDB (PDB ID 6M17, chain E, resolution 2.9 Å) (Yan et al., 2020). The advantage of the cryo-EM structures, compared to the X-ray structures, is that the fast freezing of the biological sample directly from the solution maintains the macromolecule in its closer-to-native state (Wang & Wang, 2017). The protonation state of the residues at pH = 7 was checked by using the PROPKA3.1 software (Olsson et al., 2011; Søndergaard et al., 2011), and the missing hydrogens were added by using the LEaP module in AmberTools19 (Case et al., 2019). Subsequently, the pdbqt-file of the protein target was prepared in AutoDockTools 1.5.6 (Sanner, 1999). The interactions between the RBD and the ACE2 were determined using the coordinates of ACE2 from the PDB entry 6M17 (chain B).

The 248 BLDPharm drugs that are related to the virus research () were searched by using their BD codes or CAS numbers in the ZINC15 database (Sterling & Irwin, 2015) or the PubChem database (Kim et al., 2019). Files in mol2 or sdf formats were downloaded. For some of the drugs, there were several stereoisomers in ZINC15; therefore, the total number of ligand files was increased to 327. The files in sdf format were converted to the mol2 format by using Open Babel 2.4.1 (O’Boyle et al., 2011). Finally, the pdbqt ligand files were prepared from the mol2 files by using the prepare_ligand4.py script of AutoDockTools 1.5.6 (Sanner, 1999).

Table 1.

Top-ranked ligands of the targeted docking to the RBD.

Rank	ZINC_ID	Vina energy, kcal/mol	Name	BD index	M, g/mol	logP
1	ZINC3984030	−8.5	Amentoflavone	BD113814	538.464	5.134
2	ZINC150338819	−8.4	Ledipasvir	BD290997	889.017	8.607
3	ZINC197964623	−8.3	Paritaprevir	BD304034	765.893	3.637
4	ZINC150588351	−8.2	Elbasvir	BD630189	882.035	8.116
5	ZINC85540268	−7.9	Simeprevir	BD306053	749.956	6.094
6	ZINC253632968	−7.9	Simeprevir stereoisomer 1	BD306053	749.956	5.254
7	ZINC203686879	−7.9	Velpatasvir	BD629142	883.019	7.732
8	ZINC150656835	−7.9	Simeprevir stereoisomer 2	BD306053	749.956	5.254
9	ZINC3973984	−7.9	Sotrastaurin	BD559085	438.491	2.43
10	ZINC936069565	−7.8	Glecaprevir	BD768426	838.878	3.857
11	ZINC254124762	−7.7	α-Glycyrrhizin	BD56100	822.942	2.246
12	NO ZINC ID	−7.7	Pibrentasvir	BD766533	1113.2	7.4
13	ZINC100386805	−7.6	Montelukast analog 1(same CAS number)	BD140209	572.17	8.7
14	ZINC28824700	−7.6	SC75741	BD764821	565.684	5.891
15	ZINC61389370	−7.5	Montelukast analog 2(same CAS number)	BD140209	572.17	8.7
16	ZINC68204830	−7.5	Daclatasvir	BD217334	738.89	6.222
17	ZINC220174552	−7.5	Bictegravir	BD767657	449.385	1.634
18	ZINC95551509	−7.5	Grazoprevir, MK-5172	BD293153	766.918	4.142
19	ZINC1554274	−7.4	Rilpivirine	BD211947	366.428	4.989
20	ZINC3831151	−7.4	Montelukast	BD140209	586.197	8.948

LogP – the logarithm of the octanol-water partition coefficient. This value is a measure of the hydrophobicity of a compound.

Molecular docking

The molecular docking procedure was performed in AutoDock Vina 1.1.2 (Trott & Olson, 2010). The search space was a rectangular box 40 Å × 40 Å × 28 Å that was centered at (170, 117, 240). In the PDB ID 6M17 complex, this box covers all of the ACE2 heavy atoms that are within 6 Å of the RBD (Figure 1). Due to the large size of the search space, the exhaustiveness parameter of AutoDock Vina was increased to 50.

Figure 1.

The location of the AutoDock Vina search boxes. The RBD is shown with a blue color. ACE2 is red. The atoms of ACE2 that are within 6 Å of the RBD are shown with red van der Waals spheres. The search space in the ACE2-binding region is indicated by a black rectangular box. The search space for the blind docking is indicated by a blue rectangular box.

The 20 top-ranked ligands were then subjected to the blind docking to check whether the ACE2 contact surface is their preferred binding site on RBD. The search space for the blind docking was a rectangular box 70 Å × 70 Å × 80 Å that was centered at the geometric center of the RBD (175.7, 108.4, 252.5) (Figure 1).

The visualization of the docking results was performed in VMD-1.9.3 (Humphrey et al., 1996). The 3D and 2D contact maps between the ligands and the binding site residues of the RBD were built in the PLIP tool (Salentin et al., 2015) and ligplot+ (Laskowski & Swindells, 2011), respectively. The maximum acceptable distance between the hydrogen bond (HB) acceptor (A) and the donor (D) was set at 3.35 Å. The minimum HB angle (D–H⋅⋅⋅A) was equal to 90°. The hydrophobic interactions were defined for the hydrophobic atoms (C or S) that were within 3.9 Å of each other.

Molecular dynamics simulations

To validate the docking results, the explicit solvent MD simulations for 5 top docking candidates (paritaprevir, ledipasvir, simeprevir ZINC150656835, simeprevir ZINC85540268, and amentoflavone) were performed using AMBER18 software (Case et al., 2019). For paritaprevir, an additional simulation with the neutral form of the ligand was done. Since the binding site of amentoflavone did not overlap with the binding sites of paritaprevir (anion and neutral), ledipasvir, and stereoisomers of simeprevir, the ternary complexes of RBD and two ligands (one of which was amentoflavone) were built and the MD simulations of them were also performed. A total of 12 systems were studied by MD: 6 binary ligand–RBD complexes, 5 ternary ligand–amentoflavone–RBD complexes, and the free RBD (). The partial charges of ligands were calculated using the RESP methodology in the REDIII package (Dupradeau et al., 2010). The FF14SB (Maier et al., 2015) and GAFF2 (second generation of GAFF; Wang et al., 2004, version 2.11, 2016) force fields were used for the protein and ligands, correspondingly. Each complex was put in a truncated octahedral periodic box with TIP3P (Jorgensen et al., 1983) waters (the thickness of the water shell around the solute was 15 Å). The chlorine ions were added to neutralize the system. The electrostatic interactions were treated using the particle-mesh Ewald method (Darden et al., 1993) with a direct space cutoff of 9 Å. The SHAKE method (Ryckaert et al., 1977) was applied to constrain all bonds involving hydrogen atoms. A 2 fs time integration step was used. The system was equilibrated using a multistage protocol (). The last stage of the equilibration was a 10 ns NPT MD (T = 300 K, p = 1 bar) without restraints using the Langevin thermostat (Loncharich et al., 1992) with a collision frequency of 2 ps−1 and the Berendsen barostat (Berendsen et al., 1984) with a relaxation time of 5 ps. It was followed by a 120 ns NPT MD production run with the same parameters. The snapshots were saved at every 5 ps giving a total of 24,000 frames in the MD trajectory. The analysis of the MD trajectories was done in the CPPTRAJ module of the AmberTools 19 package (Case et al., 2019). HBs were determined using the following geometric criteria: the donor–acceptor distance should be less than 3.2 Å and the donor–hydrogen–acceptor angle more than 120°. Only heteroatoms were considered as HB donors and acceptors.

Table 2.

Blind-docking-rescored 20 top-ranked ligands.

Rank	Vina energy, kcal/mol	Name
1	−9.6	Paritaprevir (anion)
2	−9.2	Ledipasvir
3	−9.1	Simeprevir stereoisomer 2
4	−8.5	Simeprevir
5	−8.5	Amentoflavone
6	−8.5	Velpatasvir
7	−8.4	Glecaprevir
8	−8.4	Elbasvir
9	−8.3	α-Glycyrrhizin
10	−8.2	SC75741
11	−8.0	Sotrastaurin
12	−7.8	Simeprevir stereoisomer 1
13	−7.8	Pibrentasvir
14	−7.7	Daclatasvir
15	−7.6	Bictegravir
16	−7.5	Montelukast
17	−7.5	Montelukast analog 2 (same CAS number)
18	−7.4	Rilpivirine
19	−7.4	Montelukast analog 1 (same CAS number)
20	−7.2	Grazoprevir, MK-5172

The ligands that were bound outside the ACE2-binding region are italicized.

Table 3.

The contacts between the ligands and the RBD in the blind-docking top-ranked complexes.

Rank	Ligand	NHB plip/ligplot+	Nπ-stacking plip	Nhydrophobic plip/ligplot+	Nbind. site res plip/ligplot+	RBD-binding site residues plip/ligplot+	RBD-binding site residues ligplot+
1	Paritaprevir	4/4	1	5/22	6/7	LEU452, LEU455, VAL483, GLU484, PHE490, GLN493	LEU452, LEU455, VAL483, GLU484, TYR489, PHE490, GLN493
2	Ledipasvir	2/3	0	8/34	8/12	ARG346, ALA348, TYR351, LEU452, GLU484, TYR489, PHE490, GLN493	ARG346, ALA348, SER349, TYR351, ALA352, ASN450, LEU452, PHE456, GLU484, TYR489, PHE490, GLN493
3	ZINC150656835 (Simeprevir stereoisomer)	2/1	2	6/32	5/6	LEU452, PHE456, TYR489, PHE490, GLN493	LEU452, PHE456, TYR489, PHE490, GLN493, SER494
4	SimeprevirZINC85540268	1/1	2	5/29	5/5	LEU452, PHE456, TYR489, PHE490, GLN493	LEU452, PHE456, TYR489, PHE490, GLN493
5	Amentoflavone	5/5	1	3/22	5/7	TYR453, GLY496, GLN498, ASN501, TYR505	ARG403, TYR453, TYR495, GLY496, GLN498, ASN501, TYR505

The RBD residues that have contact with ACE2 are marked red.

To identify the most significant structures, a cluster analysis of trajectories was performed. A hierarchical agglomerative algorithm with a 5 Å tolerance criterion was used. For the representative structure of the largest cluster, 3D and 2D contact maps were built using the PLIP tool and ligplot + software.

MMPBSA-binding free energy estimates

For the largest cluster of each MD trajectory, the binding free energy of ligands to the RBD was computed using the MMPBSA methodology (Miller et al., 2012). Energies were calculated for every 50th frame of the largest cluster trajectory (the minimum time distance between the used frames was 250 ps). The polar part of the solvation-free energy was computed using the linear Poisson-Boltzmann equation. The nonpolar part was calculated proportionally to the solvent accessible surface area. The PARSE radii (Sitkoff et al., 1994) were used and the surface tension and the offset parameters were set to 0.00542 kcal/(mol·Å2) and 0.92 kcal/mol, respectively. Genheden’s protocol (Genheden & Ryde, 2015) was used to calculate the entropy contribution: the entropy was calculated by a normal-mode analysis; the system for the minimization consisted of a ligand, as well as water molecules and protein residues within 12 Å of the ligand; all water molecules and protein residues at distances between 8 and 12 Å from the ligand were fixed in the minimization process and their contributions to the entropy were not taken into account.

Multiple independent MD simulations

To reduce the error of the MMPBSA-binding free energies, multiple independent MD simulations were performed for several systems. The snapshots from the trajectory of the largest cluster taken at a 1 ns time interval were used as starting structures for the new simulations. The initial velocities were assigned randomly according to the Maxwell distribution at a given temperature (300 K). After a short equilibration protocol (), a 4 ns NPT production phase was used to calculate the energies by the MMPBSA method: the energies were calculated at every 100 ps, for a total of 40 frames. The obtained results were averaged over all independent MD trajectories of the complex.

Table 4.

The root-mean-square deviations for the complexes of the RBD and ligands observed in the MD simulations.

Complex of RBD with ligand/ligands		RMSDbackbonedock, Å	RMSDbackboneav, Å	RMSDligdock, Å	RMSDligav, Å	RMSFlig, Å
Binary complexes	Paritaprevir (anion) ZINC197964623	4.0 ± 0.3	1.4 ± 0.4	12.1 ± 1.1	3.5 ± 2.4	4.1
	Paritaprevir (neutral)	3.7 ± 0.5	1.9 ± 0.3	3.9 ± 0.7	2.4 ± 0.8	2.5
	Ledipasvir ZINC150338819	5.0 ± 0.5	2.1 ± 0.6	11.8 ± 3.5	6.6 ± 3.1	7.4
	Simeprevir ZINC150656835	3.7 ± 0.3	1.4 ± 0.3	10.2 ± 3.0	6.4 ± 2.0	6.6
	Simeprevir ZINC85540268	2.7 ± 0.2	1.2 ± 0.3	8.8 ± 2.7	4.5 ± 2.0	4.9
	Amentoflavone ZINC3984030	3.1 ± 0.3	1.5 ± 0.3	12.8 ± 2.7	5.4 ± 3.5	6.4
Ternary complexes	Paritaprevir (anion) ZINC197964623 +	3.8 ± 0.5	2.0 ± 0.5	14.5 ± 8.9	9.0 ± 6.7	11.2
	Amentoflavone ZINC3984030			28.6 ± 12.9	23.7 ± 12.4	26.8
	Paritaprevir (neutral) +	2.9 ± 0.4	1.6 ± 0.3	5.8 ± 1.0	2.0 ± 0.7	2.0
	Amentoflavone ZINC3984030			5.6 ± 1.2	3.1 ± 0.9	3.2
	Ledipasvir ZINC150338819 +	2.8 ± 0.3	1.7 ± 0.3	9.1 ± 2.6	4.9 ± 1.2	4.7
	Amentoflavone ZINC3984030			10.9 ± 2.0	5.5 ± 2.3	5.9
	Simeprevir ZINC150656835 +	3.7 ± 0.2	1.7 ± 0.3	10.7 ± 3.8	6.0 ± 2.9	6.4
	Amentoflavone ZINC3984030			11.0 ± 3.2	8.7 ± 2.1	8.8
	Simeprevir ZINC85540268 +	3.3 ± 0.3	1.9 ± 0.3	8.3 ± 0.9	2.7 ± 0.9	2.7
	Amentoflavone ZINC3984030			10.5 ± 5.5	9.0 ± 2.0	8.9

RMSDbackbonedock – the RMSD of the RBD backbone heavy atoms with respect to the starting (PDB) structure.

RMSDbackboneav – the RMSD of the RBD backbone heavy atoms with respect to the average MD structure.

RMSDligdock – the RMSD of the ligand heavy atoms calculated with respect to the position of ligand in the docking complex, the structures of complexes were superimposed using the protein backbone heavy atoms.

RMSDligav – the RMSD of the ligand heavy atoms with respect to the position of ligand in the average MD structure, the structures of complexes were superimposed using the protein backbone heavy atoms.

RMSFig – the root-mean-square fluctuations of the ligand heavy atoms in the MD trajectory.

Results and discussion

At the first stage of the study, the targeted docking was performed to find the ligands, which block the ACE2-binding site on the RBD. The 20 top-ranked ligands are listed in Table 1. The free energies of binding for the ligands (performed via AutoDock Vina) were in a range of −8.5 to −7.4 kcal/mol. Most of the top-ranked ligands are large hydrophobic molecules. The molecular weights are greater than 700 g/mol for 8 of the first 10 top-ranked ligands, with the exception of amentoflavone (M = 538 g/mol) and sotrastaurin (M = 438 g/mol).

To check whether the ACE2-binding region on RBD is the preferred binding site for the ligands, the blind docking was performed for the 20 top-ranked ligands. The results of the blind docking are given in Table 2. The binding free energies ranged from −9.6 to −7.2 kcal/mol. For most ligands, more favorable binding free energies were obtained indicating that the small search space was not enough to find the most favorable position on the RBD. At the same time, 7 of the 20 top-ranked ligands (glecaprevir, α-glycyrrhizin, SC75741, bictegravir, montelukast, rilpivirine, grazoprevir) were bound outside of the ACE2-binding site on the RBD.

Nine of the 13 top-ranked ligands that were still bound in the ACE2 binding region are drugs used for hepatitis C virus (HCV) treatments. Paritaprevir (rank 1) and simeprevir (rank4, rank 3 [stereoisomer ZINC150656835] and rank 12 [stereoisomer ZINC253632968]) inhibit the HCV protease complex that is comprised of NSP 3 and 4A (NS3/NS4A) (McCauley & Rudd, 2016). Ledipasvir (rank 2), velpatasvir (rank 6), elbasvir (rank 8), pibrentasvir (rank 13), and daclatasvir (rank 14) are inhibitors of the NS5A protein of HCV (Gottwein et al., 2018). For all of these drugs, the successful docking results against one or several SARS-CoV-2 proteins have been previously reported in the literature (Alamri et al., 2020; Anwar et al., 2020; Balasubramaniam & Shmookler Reis, 2020; Chen, Zhang, et al., 2020; Hosseini & Amanlou, 2020; Khan, Jha, Amera, et al., 2020; Khan, Jha, Singh, et al., 2020; Maffucci & Contini, 2020; Mevada et al., 2020; Trezza et al., 2020). Several clinical trials of the combination of daclatasvir and sofosbuvir in the treatment of COVID-19 patients have been registered in Iran (Abbaspour, 2020; Foroughi, 2020; Mobarak, 2020; Roozbeh, 2020; Sadeghi, 2020). Two clinical trials using the sofosbuvir/ledipasvir combination for the treatment of COVID-19 have been completed in Egypt but no results have been posted yet (El-Gohary, 2020; Medhat et al., 2020). There has also been one clinical study in which the HCV drug danoprevir was used in combination with ritonavir for the treatment of a small group of COVID-19 patients (11 people), and all of the patients recovered (Chen, Yiu, et al., 2020). These results are encouraging for the use of HCV drugs in the treatment of COVID-19.

Amentoflavone (rank 5) is a biflavonoid found in many plants including Ginkgo biloba and Hypericum perforatum, which possesses antioxidant, antivirus, and anti-inflammatory properties and is broadly used in traditional medicine (Yu et al., 2017). It is known to be an inhibitor of CYP3A4, CYP2C9 (Kimura et al., 2010), and human cathepsin B (Pan et al., 2005). Besides, amentoflavone is an effective inhibitor (IC50 = 8.3 µM) of the SARS-CoV main protease (Ryu et al., 2010). In the recent in silico studies (Ghosh et al., 2020; Lokhande et al., 2020; Mishra et al., 2020; Patil et al., 2021; Peterson, 2020; Puttaswamy et al., 2020; Saravanan et al., 2020; Swargiary et al., 2020), amentoflavone has been reported to bind tightly to the main protease of SARS-CoV-2. There are also three docking studies targeting spike glycoprotein of SARS-CoV-2, in which amentoflavone was found among the hit compounds with the binding energies: −7.6 kcal/mol (Wei et al., 2020), −8.7 kcal/mol (Puttaswamy et al., 2020), and −10.2 kcal/mol (Rameshkumar et al., 2021). However, the binding site of amentoflavone was different in these studies. Wei and Rameshkumar obtained binding of amentoflavone outside of the ACE2-binding region. In the work of Puttaswamy et al., amentoflavone was bound in the ACE2-binding region, as well as in our study, and their result (−8.7 kcal/mol) is in good agreement with the −8.5 kcal/mol obtained for amentoflavone in this study.

To analyze the specific forces that contribute to ligand–receptor affinity, the 3D and 2D contact maps were built in PLIP tool and ligplot + for the 5 blind top-ranked ligands, which were bound in the ACE2-binding region (, S8a, and S9a). The interactions are summarized in Table 3. Since PLIP and ligplot + determine HBs and hydrophobic interactions differently, two sets of values are provided. It can be observed that ligand–RBD complexes are stabilized due to hydrophobic interactions, HBs, and π-stacking. PHE490 or TYR505 residues of the RBD are typically involved in π-stacking interactions with ligands. For most ligands, the number of hydrophobic contacts is greater or comparable to the number of HBs.

The number of the RBD residues that interacted with the ligands varied from 5 to 12, according to the ligplot + data (Table 3). To compare the binding sites of the ligands with the ACE2-binding site on the RBD (Figure 2), we built the 2D contact map of the RBD and ACE2 by using the coordinates of PDBID 6m17 (Yan et al., 2020) (Figure 3). One can observe that ACE2 interacts with seven residues of the RBD via two HBs (with RBD residues GLN493 and GLY502) and hydrophobic contacts (with RBD residues GLY476, PHE486, TYR489, GLN498 and TYR505). None of the top-ranked ligands could overlap all of the residues of the ACE2-binding site on the RBD.

Figure 2.

The interaction of the spike glycoprotein RBD (blue + yellow) with ACE2 (magenta). The residues of the RBD that are in contact with ACE2 are shown in yellow. The coordinates were taken from the PDB ID 6m17.

Figure 3.

The 2D contact map of the spike glycoprotein RBD (chain E, bottom part of the figure) and ACE2 (chain B, upper part of the figure) was built in ligplot+. The coordinates were taken from the PDB ID 6m17.

Most of the top-ranked ligands could occupy two residues of the RBD, which ACE2 binds to (Table 3). It is interesting that the HCV drugs (paritaprevir, ledipasvir, simeprevir stereoisomers) and amentoflavone have different binding sites on the RBD and block different residues of the ACE2-binding site (Table 3, Figure 4). The HCV drugs interact via HBs with GLN493 and via hydrophobic interactions with TYR489. Amentoflavone makes hydrophobic contacts with GLN498 and TYR505. Thus, their combined use is possible. Besides, when considering that amentoflavone is a CYP3A4 and CYP2C9 inhibitor, its use in combination with the HCV drugs could increase their plasma concentration.

Figure 4.

The ternary complex of the spike glycoprotein RBD (blue) with paritaprevir (red) and amentoflavone (yellow).

In agreement with our results, two different binding sites for ligands were identified on RBD in the ACE2 contact region in several other studies (Hakmi et al., 2020; Maffucci & Contini, 2020). The first binding site is composed of LEU455, PHE456, PHY486, ASN487, TYR489 and GLY493 RBD residues and corresponds to the HCV drug-binding site in our case. The second site (formed by Tyr449, GLY496, GLN498, THR500, ASN501, GLY502 and Tyr505 residues) is the binding region of amentoflavone.

Using the blind docking poses of the 5 top-ranked ligands we built the ternary complexes of the RBD, amentoflavone, and the HCV drugs (Figure 4, ). For paritaprevir, the neutral form of the ligand was also considered. The blind docking pose of the neutral paritaprevir was similar to the anion pose (Supplementary Figures S1b,c). The result of the blind docking for the neutral form of paritaprevir was −9.2 kcal/mol.

To check the stability of the docking complexes, the MD simulations of the binary and ternary ligands–RBD complexes were performed. Table 4 shows the root-mean-square deviations (RMSD) for the MD trajectories of the complexes calculated relative to the starting (docking) structures as well as to the average MD structures. The RMSDs of the RBD backbone with respect to the average structures did not exceed 2.1 Å, which indicates that the protein structure was relatively stable during the MD simulations. At the same time, the RMSDs of the RBD backbone relative to the starting structure of the protein (PDBID 6m17) ranged from 2.7 to 5.0 Å for the trajectories of different complexes (3.3 Å for the RBD trajectory without ligands). It is a sign of protein conformational changes during the MD simulation. The analysis of the root-mean-square fluctuations (RMSF) of the RBD residues showed that the most mobile was the RBD loop formed by residues 474–487 (Supplementary Figure S10).

The RMSDs of ligands calculated with respect to their positions in the starting docking complexes were greater than 8 Å for all complexes, except for the complexes of the neutral form of paritaprevir. This means that the structure of the docking complexes changed significantly during the MD simulation. However, in most cases, ligands remained bound to the RBD. The dissociation of ligands was observed only in the MD trajectory of the ternary complex of the anionic form of paritaprevir, amentoflavone and RBD. The HCV drugs (the neutral form of paritaprevir, ledipasvir, and both stereoisomers of simeprevir) had lower RMSF values in the ternary complexes compared to the binary complexes, meaning that they were less mobile in the ternary complexes. According to the RMSD and RMSF data, the complexes of the neutral form of paritaprevir were the most stable: they were characterized by the lowest RMSF (2.5 Å and 2.0 Å) and RMSDligdock (3.9 Å and 5.8 Å) values (Table 4).

Due to the relatively flat RBD–ACE2-binding surface, it is challenging to find ligands that will remain tightly bound to it during the MD simulation. In several other MD studies of ligand–RBD complexes, the dissociation or migration of ligands to the other binding sites was observed in the course of the simulation (Chikhale et al., 2020; Deganutti et al., 2021; Maffucci & Contini, 2020; Unni et al., 2020).

Analysis of contacts between the ligands and the RBD in the MD trajectories

Given the mobility of ligands in complexes during the MD simulation, a cluster analysis of trajectories was performed to identify the most significant structures. The representative structures of the main MD cluster for each complex are shown in Supplementary Figures S2–S6b,c and S7. For the representative structures, the analysis of contacts between the ligand and the RBD was performed (Table 5, Supplementary Figures S2–S6e,f,h,i, S8, and S9). The stabilization of complexes was due to the HBs and hydrophobic interactions. In the binary ligand–RBD complexes, most ligands (the neutral form of paritaprevir, ledipasvir, both simeprevir stereoisomers) remained bound in the region of the ACE2-binding site during the MD simulation and interacted with the RBD residues GLN493 and TYR489. At the same time, amentoflavone and the anion form of paritaprevir moved along the RBD to the new binding sites that did not overlap with the ACE2-binding site. Interestingly, in most ternary complexes, amentoflavone, due to the interaction with the second ligand, remained within the ACE2-binding region during MD and blocked TYR505 and GLN498 residues of the RBD (which bind to ACE2). This gives evidence in favor of the combined use of the HCV drugs and amentoflavone to block the binding of ACE2 to RBD.

Table 5.

The contacts between the ligands and the RBD in the representative structures of the main clusters of MD trajectories.

Complex of RBD with ligand/ligands		NHB plip/ligplot+	Nπ-stack. plip	Nhydroph. plip/ligplot+	Nbind. site res plip/ligplot+	RBD-binding site residues plip	RBD-binding site residues ligplot+
Binary complexes	Paritaprevir (anion)	0/0	0	5/16	2/7	TYR449, PHE490	GLY446, GLY447, TYR449, ASN450, LEU452, PHE490, SER494
	Paritaprevir (neutral)	2/2	0	3/17	3/5	TYR489, PHE490, GLN493	ILE472, TYR489, PHE490, GLN493, SER494
	Ledipasvir	4/8	1	5/29	7/10	TYR449, LEU452, GLY485, TYR489, PHE490, GLN493, SER494,	TYR449, LEU452, PRO479, GLY485, ASN487, CYS488, TYR489, PHE490, GLN493, SER494
	Simeprevir ZINC150656835	2/3	1	ARG403, LEU455, PHE456, TYR473, TYR489, GLN493	ARG403, TYR453, LEU455, PHE456, TYR473, TYR489, GLN493, TYR495	ARG403, TYR453, LEU455, PHE456, TYR473, TYR489, GLN493, TYR495
	Simeprevir ZINC85540268	3/3	0	4/14	6/7	TYR351, TYR449, LEU452, LEU492, GLN493, SER494	TYR351, TYR449, ASN450, LEU452, LEU492, GLN493, SER494
	Amentoflavone	2/1	0	3/25	3/6	GLY404, VAL503, TYR508	GLY404, ASP405, ARG408, VAL503, GLY504, TYR508
Ternary complexes	Paritaprevir (anion) +	0/0	0	4/10	3/3	TYR449, THR470, PHE490	TYR449, THR470, PHE490
	Amentoflavone	–	–	–	–	–	–
	Paritaprevir (neutral) +	2/2	0	7/30	5/8	LEU452, ILE472, TYR489, PHE490, LEU492	LEU452, THR470, ILE472, GLU484, GLY485, TYR489, PHE490, LEU492
	Amentoflavone	2/3	0	1/10	3/5	TYR495, ASN501, TYR505	TYR495, GLY496, GLN498, ASN501, TYR505
	Ledipasvir +	4/5	0	4/19	6/9	TYR449, LEU452, GLU471, PRO491, LEU492, GLN493	TYR449, LEU452, THR470, GLU471, PRO491, LEU492, GLN493, SER494, UWI
	Amentoflavone	3/2	0	0/3	2/3	TYR449, ASN501	TYR449, ASN501, RTW
	Simeprevir ZINC150656835 +	1/1	1	4/39	3/6	ILE472, PHE490, GLN493	TYR449, ILE472, PHE490, GLN493, SER494, UWI
	Amentoflavone	2/2	0	1/16	2/4	TYR449, TYR505	TYR449, ASN501, TYR505, SQY
	Simeprevir ZINC85540268 +	0/0	0	3/16	2/6	VAL483, PHE490	TYR449, GLY482, VAL483, CYS488, PHE490, LEU492
	Amentoflavone	1/4	1	1/7	3/6	GLN498, ASN501, TYR505	SER494, TYR495, GLY496, GLN498, ASN501, TYR505

NHB plip/ligplot+ – the number of hydrogen bonds between the ligand and the RBD in the plip/ligplot + contact maps.

Nπ-stack. plip – the number of π-stacking interactions between the ligand and the RBD in the plip contact map.

Nhydroph. plip/ligplot+ – the number of hydrophobic interactions between the ligand and the RBD in the plip/ligplot + contact maps.

Nbind. site res plip/ligplot+ – the number of the RBD residues that interact with the ligand according to the plip/ligplot + contact maps.

The RBD residues that have contact with ACE2 are bold.

The dynamics of the intermolecular HBs in the ligands–RBD complexes was analyzed using the MD trajectories (). The average number of HB between the ligand and the RBD during the MD simulation ranged from 0 to 3. The highest average number of the intermolecular HB was observed for ledipasvir (2.8) in the binary complex and the neutral form of paritaprevir in the ternary complex (3.0) (Supplementary Figures S11 and S12, ). Since the binding sites of the HCV drugs and amentoflavone were different, they formed HB with the different residues of the RBD. The functional groups of the HCV drugs that participated in hydrogen bonding with the RBD during the MD simulation included: carboxamide groups (ledipasvir, paritaprevir, simeprevir), sulfonamide group (simeprevir, paritaprevir), benzimidazole and imidazole rings (ledipasvir), methoxy group (ledipasvir), and fluorine (ledipasvir). Most ligands (HCV drugs) formed HB with the GLN493 residue of the RBD, which interacts via HB with ACE2 (Figure 3). The other RBD residues that were in contact via HB with the HCV drugs were: ARG403, TYR453, GLU471, GLY485, ASN487, CYS488, PHE490, LEU492, and SER494.

Amentoflavone formed HB with the RBD via hydroxyl groups. In the trajectory of its binary complex, the HB contacts with the RBD residues GLY404 and GLY504 were detected. In the ternary complexes, amentoflavone participated in hydrogen bonding most frequently with the ASN501 residue of RBD, among other residues were ASP405, TYR449, SER494, TYR495, GLY496, GLN498, GLY502. Many of these residues are important for the RBD–ACE2 interaction (Yan et al., 2020).

The average lifetime of the intermolecular HB in the ligands–RBD complexes was 10–100 ps.

To estimate the binding affinities of ligands to the RBD, the MMPBSA-binding free energies were calculated for the main MD clusters of the trajectories of complexes (Table 6). In the binary ligand–RBD complexes, the negative (favorable) binding free energy was obtained only for ledipasvir (−5.9 kcal/mol). The complexes were stabilized by van der Waals, electrostatic and nonpolar solvation contributions, whereas polar solvation and entropy terms were unfavorable. The analysis of the entropy components () showed that the main unfavorable contributions come from the loss of translational and rotational degrees of freedom of the ligand. For all the ligands, this penalty is circa 25 kcal/mol. The entropy term seems to be the most challenging in the MMPBSA calculation and many studies ignore it. At the same time, the omitting of the entropy contribution could worsen the results (Genheden et al., 2012). In this study, the translational and rotational entropy terms were calculated using standard statistical mechanics formulas for rigid rotor-harmonic oscillator-ideal gas approximation (McQuarrie, 1976). This approach is known to overestimate the translational entropy loss in the case of the liquid phase (Mammen et al., 1998). So, the calculated entropy penalty is probably overestimated, which resulted in the positive binding free energies for several ligands. Still, we suppose that the relative ranking of ligands was not affected. It should be noted that all the ligands had favorable enthalpies of binding, confirming their possible interaction with RBD.

Table 6.

The MMPBSA-binding free energies for the ligand–RBD complexes.

Ligand	Nfr 1-st cl/% of traj time	Nframes MMPBSA	VdW, kcal/mol	EEL, kcal/mol	EPB, kcal/mol	ENPOLAR (SASA), kcal/mol	dH, kcal/mol	−TdS, kcal/mol	dG, kcal/mol
Paritaprevir (anion)ZINC197964623	19,688/82.0	394	−37.5 ± 4.1	−36.9 ± 7.3	57.6 ± 7.4	−4.7 ± 0.4	−21.4 ± 3.2	26.5 ± 5.0	5.1 ± 5.9
Paritaprevir(neutral)	23,200/96.7	464	−36.2 ± 3.9	−13.3 ± 5.7	37.4 ± 7.3	−5.0 ± 0.4	−17.1 ± 3.5	23.8 ± 4.4	6.7 ± 5.6
LedipasvirZINC150338819	10,532/43.9	211	−58.4 ± 6.0	−31.4 ± 10.0	63.5 ± 9.6	−6.6 ± 0.4	−32.9 ± 6.0	27.0 ± 5.0	−5.9 ± 7.8
SimeprevirZINC150656835	14,291/59.5	286	−40.3 ± 4.8	−110.2 ± 15.1	133.7 ± 15.2	−5.4 ± 0.4	−22.2 ± 4.4	25.9 ± 5.0	3.7 ± 6.7
SimeprevirZINC85540268	8623/35.9	173	−39.5 ± 3.1	−64.1 ± 11.8	83.8 ± 11.7	−5.1 ± 0.4	−24.9 ± 3.3	27.2 ± 4.6	2.3 ± 5.7
AmentoflavoneZINC3984030	17,494/72.9	350	−29.6 ± 2.5	−14.3 ± 4.2	29.4 ± 3.6	−4.0 ± 0.1	−18.4 ± 2.9	18.8 ± 3.0	0.4 ± 4.2
Paritaprevir anion + (Amentoflavone)a	6806/28.4	137	−35.6 ± 5.9	−35.9 ± 9.8	58.1 ± 13.0	−4.8 ± 0.5	−18.2 ± 3.1	24.9 ± 4.6	6.7 ± 5.5
Paritaprevir neutral +	24,000/100.0	480	−55.4 ± 4.7	−27.9 ± 8.1	57.7 ± 9.5	−6.6 ± 0.4	−32.3 ± 3.6	23.9 ± 5.3	−8.4 ± 6.4
Amentoflavone	Amentoflavone		−29.2 ± 3.9	−16.1 ± 6.7	35.8 ± 7.7	−4.1 ± 0.4	−13.6 ± 3.7	20.4 ± 4.1	6.8 ± 5.5
Ledipasvir +	10,928/45.5	219	−55.5 ± 5.7	−16.9 ± 4.8	41.8 ± 4.7	−6.3 ± 0.4	−36.9 ± 5.1	27.8 ± 4.4	−9.1 ± 6.7
	Amentoflavone		−27.3 ± 3.6	−15.4 ± 9.8	34.2 ± 10.1	−4.1 ± 0.3	−12.6 ± 3.4	18.6 ± 3.8	6.0 ± 5.1
Simeprevir ZINC150656835 +	8761/36.5	176	−38.8 ± 5.1	−55.8 ± 10.2	72.8 ± 11.9	−5.0 ± 0.5	−26.8 ± 3.4	25.2 ± 4.7	−1.6 ± 5.8
	Amentoflavone		−39.1 ± 4.9	−21.5 ± 5.3	50.6 ± 6.1	−5.1 ± 0.3	−15.1 ± 3.8	19.3 ± 4.8	4.2 ± 6.1
SimeprevirZINC85540268 +	12,808/53.4	257	−40.7 ± 3.7	−46.8 ± 8.8	65.2 ± 9.0	−5.0 ± 0.3	−27.4 ± 3.4	24.5 ± 5.2	−2.9 ± 6.2
	Amentoflavone		−27.0 ± 3.8	−27.7 ± 10.2	46.9 ± 10.2	−4.1 ± 0.3	−11.8 ± 4.5	22.1 ± 4.5	10.3 ± 6.4

aThere are no data for amentoflavone because it dissociated from this ternary complex during the MD simulation.

Nfr 1-st cl – the number of frames in the first (largest) cluster.

Nframes MMPBSA – the number of frames that were used for the MMPBSA calculation.

VdW – the van der Waals contribution to the binding free energy.

EEL – the electrostatic contribution to the binding free energy.

EPB – the polar contribution to the solvation-free energy.

ENPOLAR – the nonpolar contribution to the solvation free energy; ENPOLAR = γ·SASA + b, where SASA – is the solvent-accessible surface area, γ = 0.00542 kcal/(mol·Å2) and b = 0.92 kcal/mol.

dH – the enthalpy contribution to the binding free energy; dH = VdW + EEL + EPB + ENPOLAR.

−TdS – the entropy contribution to the binding free energy.

dG – the binding free energy; dG = dH – TdS.

In the ternary complexes, all ligands had the negative binding free energies except for the anion form of paritaprevir and amentoflavone (which was present in each ternary complex as the second ligand). The affinity of ligands to the RBD in the ternary complexes decreases in a row: ledipasvir (−9.1 kcal/mol) > neutral form of paritaprevir (−8.4 kcal/mol) > simeprevir ZINC85540268 (−2.9 kcal/mol) > simeprevir ZINC150656835 (−1.6 kcal/mol). Based on the MMPBSA energy estimates, it can be concluded that ledipasvir is the best candidate for binding to the RBD. Although amentoflavone displays weak binding to the RBD, its presence promotes the stronger binding of other ligands (HCV drugs) in the ternary complexes.

To increase the precision of the MMPBSA-binding free energies, additional independent MD simulations were performed for the ledipasvir–RBD, ledipasvir–amentoflavone–RBD, and neutral paritaprevir–amentoflavone–RBD complexes. The standard deviations of the binding free energies were reduced by 1.5–2 times (Table 7). The relative affinities of ligands for the RBD remained unchanged: in the ligand–amentoflavone–RBD ternary complexes, the binding free energy of ledipasvir (−9.5 kcal/mol) was more favorable than the binding free energy of the neutral form of paritaprevir (−5.8 kcal/mol). A comparison of the binary and ternary complexes of ledipasvir shows that the presence of amentoflavone leads to the formation of a stronger complex of ledipasvir with RBD: the binding free energy of ledipasvir is 2 times larger in the ternary complex.

Table 7.

The MMPBSA-binding free energies calculated from the multiple independent simulations.

Ligand		Number of simulations	dH, kcal/mol	−TdS,kcal/mol	dG, kcal/mol
Binary complex	Ledipasvir	56	−32.6 ± 3.7	27.9 ± 1.5	−4.7 ± 4.0
Ternary complexes	Ledipasvir +	55	−37.0 ± 4.3	27.5 ± 1.1	−9.5 ± 4.4
	Amentoflavone		−12.4 ± 1.4	18.9 ± 1.0	6.5 ± 1.7
	Paritaprevir neutral +	120	−31.6 ± 2.2	25.8 ± 1.3	−5.8 ± 2.6
	Amentoflavone		−12.9 ± 2.5	19.5 ± 1.3	6.6 ± 2.8

The comparison of the MMPBSA-binding energies with the literature data is difficult because the results strongly depend on the parameters and details of the simulation (Wang et al., 2019). The choice of the polar solvation energy calculation method has a significant impact on the results: the MMGBSA-binding free energies are too negative compared to the MMPBSA ones. Due to the fact that in many studies the entropy contribution is neglected, we compared the enthalpy, which was calculated as the sum of the van der Waals, electrostatic, polar (PB), and nonpolar (SASA) solvation energy terms, with the binding energies reported in the literature. The enthalpies for the ligand–RBD binary complexes obtained in this study varied from −32.9 to −17.1 kcal/mol (Table 6). Our results compare well with the reported −18.4 kcal/mol for paromycin–RBD complex (Tariq et al., 2020), −19.6 kcal/mol for apigenin-7-O-rutinoside (Albohy et al., 2020), −17.9 kcal/mol for vancomycin and −20.9 kcal/mol for teicoplanin (Tao et al., 2021), and −15.6 kcal/mol for digitoxin (Wei et al., 2020). de Oliveira et al. (2020) obtained slightly more favorable energies for theaflavin digallate and suramin: −38.5 kcal/mol and −40.4 kcal/mol, respectively. In the study of Trezza et al. (2020), a stable MD trajectory was obtained for simeprevir–RBD complex. They calculated the simeprevir–RBD interaction energy as a sum of short-range electrostatic and van der Waals terms and got −18 kcal/mol.

The reported in the literature MMGBSA energies for ligand–RBD complexes are, in general, more negative: −89.4 kcal/mol for withanoside X (Chikhale et al., 2020), −52.1 kcal/mol for the top-ranked ligand 5960 in the study of Hakmi et al. (2020); −53.2 kcal/mol for nilotinib (Deganutti et al., 2021). Maffucci used 60 explicit water molecules during the calculation of the MMGBSA energies for ligand–RBD complexes (Maffucci & Contini, 2020). They obtained −164.3 kcal/mol for the top-ranked ligand polymyxin B and −60.8 kcal/mol for ledipasvir, which was among the top 2% of compounds.

Thus, the literature data confirm our results on the stability of the ledipasvir–RBD and simeprevir–RBD complexes. To the best of our knowledge, there are no other MD simulations of paritaprevir–RBD and amentoflavone–RBD complexes. But, there are many stable MD simulations of the amentoflavone complex with the main protease (Ghosh et al., 2020; Lokhande et al., 2020; Mishra et al., 2020; Patil et al., 2021; Saravanan et al., 2020; Swargiary et al., 2020). For paritaprevir and simeprevir, stable MD simulations with the main protease and the uridylate-specific endoribonuclease (NSP15) have been obtained (Alamri et al., 2020; Khan, Jha, Amera, et al., 2020; Khan, Jha, Singh, et al., 2020). Ledipasvir has been also shown to bind to the main protease, in addition to the spike RBD (Maffucci & Contini, 2020). So, for all of our top-ranked ligands, the activity against other SARS-CoV-2 protein targets has been reported in the literature.

Summarizing the MMPBSA results of this study, ledipasvir and the neutral form of paritaprevir when used in combination with amentoflavone are the most efficient RBD binders. The representative structures of the main MD clusters of the ternary complexes of ledipasvir and the neutral form of paritaprevir are shown in Figures 5 and 6, correspondingly.

Figure 5.

The representative structure of the main MD cluster of the ledipasvir–amentoflavone–RBD complex.

Figure 6.

The representative structure of the main MD cluster of the neutral paritaprevir–amentoflavone–RBD complex.

Conclusions

Due to the fast and uncontrolled spread of the COVID-19 disease caused by the SARS-CoV-2 virus, there is a need to find a safe and effective treatment. Since the virus penetrates in the host cell with the help of spike glycoprotein, this protein represents an attractive target for drug repurposing. We performed the molecular docking of 248 drugs against the RBD of the spike glycoprotein of SARS-CoV-2. Most of the top-ranked ligands are the FDA-approved HCV drugs: paritaprevir, ledipasvir, simeprevir. In the top 5, there is also a natural product, the biflavonoid amentoflavone. The HCV drugs and amentoflavone have different binding sites on the RBD, thus, the ternary complexes of two ligands and RBD can be created.

To validate the docking results, the binary and ternary ligands–RBD complexes were studied by MD simulation with explicit solvent for 120 ns. The structure of the docking complexes changed during the simulation, but in most cases, ligands remained within the ACE2-binding region. The smallest structural changes relative to the starting docking configuration were observed for the complexes of the neutral form of paritaprevir.

According to the MMPBSA-binding free energy estimates, ledipasvir (−4.7 ± 4.0 kcal/mol) has the best result in the binary ligand–RBD complexes. The HCV drugs display stronger binding in the ternary complexes (HCV drug–amentoflavone–RBD complexes) compared to the binary complexes. The most stable ternary complexes with RBD are formed by ledipasvir (−9.5 ± 4.4 kcal/mol) and the neutral form of paritaprevir (−5.8 ± 2.6 kcal/mol).

The ligands–RBD complexes are stabilized by HBs and hydrophobic interactions. The detailed analysis of the intermolecular HBs in complexes during MD simulations showed that most HCV drugs form HBs with the RBD residue GLN493, which interacts via hydrogen bonding with ACE2. In the ternary complexes, amentoflavone remains in the ACE2-binding region during the MD simulation and forms a HB with the ASN501 residue of RBD. The MD simulations of ternary complexes confirm that simultaneous binding of the HCV drug and amentoflavone to RBD in the ACE2-binding region is possible, but experimental validation is required.

According to the literature data, the top-ranked ligands defined in this study display also activity against other proteins of SARS-CoV-2 (main protease, NSP15).

The results of our study suggest that the HCV drugs ledipasvir and paritaprevir (in a neutral form) in combination with a natural biflavonoid amentoflavone are the most promising candidates to block the binding site of ACE2 on RBD. However, one must keep in mind that the presence of the ACE2 could influence the affinity of ligands to the RBD. The next steps should be the simulation of the interaction between the ligands–RBD ternary complex and ACE2 and the validation of the results through the experimental investigations.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.